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ABSTRACT 


Four  topics  are  discussed  in  this  progress  report. 

The  first  topic  is  related  to  the  masking  of  underground  nuclear 
tests  with  large  earthquakes.  Computer  simulations,  using  sig¬ 
nals  from  actual  seismic  events,  suggest  that  test  detonations 
one  magnitude  unit  smaller  than  the  earthquake  can  be  detected, 
providing  the  shot  time  is  some  two  minutes  after  the  time  of 
occurrence  of  the  earthquake.  The  second  topics  deals  with  the 
effects  of  signal-to-noise  ratio  and  local  travel  time-anomalies 
on  the  accuracy  of  epicenter  location  using  large  seismic  arrays. 
Both  theoretical  analyses  and  computer  simulations  indicate 
standard  deviations  of  angular  errors  on  the  order  of  0.3°  for 
either  a  10  dB  signal-to-noise  ratio  or  a  peak  anomaly  on  the 
order  of  l/20th  of  a  second.  The  third  topic  deals  with  the 
application  of  a  previously  developed  coda-correlation  discriminant 
to  DIMUS  (hardlimited)  seismograms.  This  discriminant  makes  use 
of  the  average  paired  correlation  coefficient  of  the  10-second 
portion  of  seismic  arrivals  commencing  3  seconds  after  P-wave 
onset  for  an  array  of  widely  separated  stations.  Results  to 
date  indicate  a  moderate  degradation  of  the  discriminant  for 
the  DIMUS  seismograms  compared  to  the  unclipped  seismograms. 

The  fourth  topic  deals  with  the  automatic  identification  of  the 
pP  phase  of  earthquakes,  for  events  in  the  40  to  150  km  depth 
region.  Excellent  results  have  been  obtained  with  this  automatic 
scheme,  and  not  only  is  the  pP  phase  properly  identified,  but 
also  the  sP  phase.  This  combination  allows  depths  to  be  estimated 
with  considerable  confidence. 
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SECTION  I 


INTRODUCTION 

This  is  the  Second  Quarterly  Technical  Report  on 
Contract  F19628-67-C-0370 .  The  report  summarizes  progress  on 
four  topics:  masking  of  underground  tests  with  large  earth¬ 
quakes,  beamforming  and  other  techniques  for  epicenter  location, 
a  DIMUS  (hardlimited)  version  of  the  coda-correlation  discrimi¬ 
nant,  and  an  extension  of  the  automatic  pP  test  to  depths  greater 
than  40  km. 

Section  II  is  a  continuation  of  a  discussion  in  our 
First  Quarterly  Technical  Report  on  the  problem  of  detecting 
covert  nuclear  tests  in  the  vicinity  of  large  earthquakes.  In 
this  discussion,  a  detection  system  consisting  of  a  network  of 
single  stations  surrounding  the  earthquake  epicenter  plus  a 
continental-size  array  was  proposed.  Under  the  assumption  that 
the  single-station  network  would  detect  a  covert  test  if  the 
signal  from  the  test  preceded  that  from  the  earthquake  at  any 
of  the  stations,  minimum  earthquake-test  delay  times  necessary 
for  the  test  to  avoid  detection  were  derived  as  a  function  of 
test  site  position.  The  continental-size  array  imposes  maximum 
earthquake-test  delay  times  for  the  test  to  avoid  detection. 

Some  calculations  were  presented  indicating  the  general  properties 
of  the  earthquake  contribution  to  the  array  output. 

In  order  to  obtain  a  measure  of  the  sensitivity  of  the 
continental-size  array,  we  have  performed  several  additional 
calculations  that  are  summarized  in  this  report.  The  most  impor¬ 
tant  of  these  are  based  on  a  simulated  detection  problem  with 
an  assumed  nuclear  test  near  a  larger  earthquake.  This  simulation 
is  based  on  actual  seismic  data,  and  involves  several  values  of 
the  following  parameters:  earthquake-test  delay  time,  test  site 
location,  relative  magnitudes  of  the  two  events.  Results  are 
summarized  in  terms  of  plots  that  give  the  maximum  delay  time 
imposed  by  the  array  as  a  function  of  test  site  location  for 
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each  of  two  test  magnitudes.  These  plots  suggest  that  if  the 
enormous  number  of  computations  necessary  to  implement  this 
detection  system  could  be  managed,  it  would  very  significantly 
reduce  the  probability  that  a  covert  nuclear  test  would  go 
undetected . 

Section  III  presents  several  theoretical  derivations 
related  to  the  problem  of  epicenter  location  by  beamforming 
or  by  plane-wave  fitting  following  individual  arrival  time 
measurements  at  each  station.  A  highly  idealized  model  of  the 
problem  is  analyzed  with  the  objective  of  gaining  insight  into 
the  actual  seismic  problem.  More  specifically,  the  plane-wave 
fit  based  on  arrival-time  measurements  is  stated,  and  three 
methods  of  arrival-time  measurement  are  analyzed.  In  general, 
the  best  of  the  arrival-time  measurements  is  the  correlation 
time  pick, which  requires  perfect  knowledge  of  the  signal  wave¬ 
shape.  The  beamforming  calculation  is  also  analyzed,  and  it 
is  shown  that  it  yields  the  same  results  as  the  correlation 
time  pick  followed  by  a  plane-wave  fit.  This  indicates  one 
important  advantage  of  the  beamforming,  since  this  calculation 
does  not  require  prior  knowledge  of  the  exact  signal  waveshapes. 
The  results  of  this  section  are  summarized  in  terms  of  a  few 
examples.  In  one  of  these, a  uniform  21-element  array  of  total 
aperture  200  km,  with  a  phase  velocity  of  20  km/sec  and  a 
signal-to-noise  ratio  of  10  dB,  leads  to  a  standard  deviation 
in  angular  measurement  of  0.3° • 

The  principal  limitation  of  the  theoretical  analysis 
is  that  closed  form  solutions  are  available  only  with  rather 
idealized  assumptions.  In  order  to  study  more  realistic  assump¬ 
tions  under  controlled  conditions,  we  have  also  performed  several 
simulations.  These  simulations,  which  are  summarized  in  Sec¬ 
tion  IV,  are  based  on  actual  seismic  records  from  a  large  surface- 
focus  event  recorded  at  LASA-Montana .  In  these  simulations 
controlled  amounts  of  additive  noise  or  travel-time  errors  are 
added  to  each  seismic  record  and  four  different  methods  of 
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epicenter  estimation  are  tested.  Two  of  these  methods  are  based 
on  conventional  beamforming:  in  one  case  the  pattern  is  defined 
by  an  energy  calculation,  in  the  other  by  a  peak  value.  Two 
calculations  based  on  DIMUS  (hardlimited)  seismograms  are  also 
presented.  In  all  cases  an  automatic  event  detector  is  used 
to  define  the  interval  of  interest  so  that  the  calculations  will 
be  as  realistic  as  possible. 

The  results  of  these  simulations  are  summarized  in 
tables  and  figures,  which  show  that  the  additive  noise  leads 
to  errors  comparable  to  those  predicted  by  the  simplified  theory 
of  Section  III,  whereas  the  simulated  time  errors  lead  to  more 
serious  errors  than  expected.  The  reason  for  the  latter  discre¬ 
pancy  is  not  well  understood,  but  is  probably  a  consequence  of 
modeling  sampled  data  with  continuous  waveforms.  There 
is  considerable  spread  in  the  simulated  data  now  available,  but 
it  appears  that  both  of  the  DIMUS-based  calculations  yield  per¬ 
formance  levels  quite  comparable  to  those  obtained  from  analog 
waveforms . 

Section  V  is  a  brief  report  on  calculations  of  the 
DIMUS  (hardlimited)  version  of  the  coda-correlation  discriminant, 
which  worked  well  with  analog  waveforms.  The  calculations  sum¬ 
marized  in  Section  V  indicate  that  the  discriminant  based  on  DIMUS 
waveforms  leads  to  a  less  complete  (but  possibly  useful)  separa¬ 
tion  of  moderate  depth  earthquakes  from  other  events. 

Section  VI  reports  experimental  results  of  applying 
an  extended  automatic  pP  test  to  deep  earthquakes,  using  data 
from  continental-size  arrays.  The  result  of  processing  eight 
events  are  very  encouraging.  In  many  cases  the  test  has  iden¬ 
tified  both  pP  and  sP.  In  some  cases  the  test  has  yielded  depth 
estimates  that  conflict  with  C&GS  depths.  As  is  illustrated 
in  Section  VI,  it  appears  in  some  cases  that  the  C&GS  depths 
result  from  mistaking  sP  for  pP .  The  test  results  are  summarized 
in  tables  and  figures,  and  several  seismograms  are  presented 
to  confirm  and  explain  the  test  results. 
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SECTION  II 


DETECTION  OF  UNDERGROUND  NUCLEAR  EXPLOSIONS  IN 
THE  PRESENCE  OF  LARGE  NATURAL  EVENTS 

The  problem  of  detecting  an  underground  nuclear  explo¬ 
sion  that  occurs  shortly  after  an  earthquake  and  in  the  same 
geographical  vicinity  was  discussed  in  an  earlier  report  [1]. 

A  hypothetical  detection  system  consisting  of  a  continental- 
size  array  plus  a  set  of  single  stations  encircling  the  earth¬ 
quake  epicenter  at  close  ranges  (perhaps  ^0°)  was  considered. 

It  was  concluded  that,  from  the  point  of  view  of  trying  to  hide 
a  nuclear  test  in  the  earthquake's  signal,  the  continental-size 
array  imposes  a  maximum  allowable  delay  between  the  earthquake 
and  the  test,  which  is  a  function  of  both  the  relative  locations 
and  magnitudes  of  the  test  and  the  earthquake.  Roughly  speaking, 
the  maximum  delay  is  the  time  after  which  the  earthquake's  signals 
have  decayed  sufficiently  that  an  underground  nuclear  test  could 
be  detected.  The  single-station  network,  on  the  other  hand, 
imposes  a  limit  on  how  soon  after  the  earthquake  a  nuclear  test 
can  be  concealed  without  being  detected  by  one  of  these  single¬ 
seismometer  stations.  This  minimum  delay  is  imposed  by  the 
assumption  that  if  the  signal  from  the  test  precedes  that  of 
the  earthquake  at  any  single  station,  the  test  will  be  detected. 

Plots  of  the  energy  received  in  nonoverlapping  one- 
second  time  windows  as  a  function  of  time-after-onset  were  pre¬ 
sented  for  several  stations  that  recorded  two  large  natural 
events.  It  was  observed  that  the  energy-decay  rate  varied  from 
station  to  station  and  that  the  energy  decay  curve  for  a  single 
station  was  far  from  smooth.  The  arrival  of  additional  phases 
after  the  P-phase,  with  varying  relative  strengths,  produces 
several  irregular  (and  often  large)  peaks  in  the  energy  vs.  time 
curve,  both  for  single  stations  and  for  the  array  output.  The 
apparently  irregular  behavior  of  these  curves  makes  it  unlikely 


that  any  simple  model  would  suffice  to  predict  detection  perfor¬ 
mance  of  a  continental-size  array  in  this  context.  For  this 
reason  we  have  resorted  to  a  simulation,  using  actual  seismic 
records  from  a  large  earthquake,  to  study  the  sensitivity  of  a 
continental-size  array  as  a  function  of  the  relative  magnitudes 
and  locations  of  the  two  events  and  the  time  between  them. 

Two  methods  of  detection,  which  were  motivated  by  the 
time-varying  nature  of  the  coda  energy  levels,  have  been  studied 
by  means  of  simulation.  These  two  methods  are  very  similar  to 
calculations  used  in  an  automatic  pP  enhancement  method  which 
was  discussed  in  an  earlier  report  [2]..  All  of  the  results 
presented  in  this  section  are  based  on  a  magnitude-5.5  earthquake 
as  recorded  at  fifteen  stations  located  in  the  continental  United 
States.  The  earthquake  occurred  in  the  Fox  Islands  on  1  Sept.  1964. 

The  contents  of  this  section  may  be  summarized  as 
follows.  The  energy-ratio  detection  statistic,  e,  which  is 
defined  as  the  average  power  in  a  one-second  window  divided  by 
that  of  the  preceding  five  seconds,  is  first  discussed.  Then 
another  detection  statistic,  which  is  obtained  by  multiplying  e 
by  the  average  paired  correlation  coefficient  obtained  for  the 
same  one-second  interval,  is  considered.  Following  the  presen¬ 
tations  of  these  two  detection  statistics,  the  expected  perfor¬ 
mance  of  both  is  studied,  and  conclusions  are  drawn  indicating 
how  the  minimum  detectable  magnitude  of  an  event  may  be  related 
to  the  average  power  in  the  output  of  the  continental-size  array. 
Following  this  discussion,  curves  of  minimum  detectable  nuclear- 
test  magnitude  as  a  function  of  the  earthquake-shot  delay  time, 
for  selected  points  in  the  vicinity  of  the  earthquake's  epicenter 
are  presented  and  discussed.  Finally,  the  results  of  this  simu¬ 
lation  are  summarized  in  terms  of  allowable  intervals  for  escaping 
detection  with  the  methods  discussed,  as  a  function  of  the  nuclear 
test's  magnitude  and  the  relative  position  of  the  test  site  and 
the  earthquake's  epicenter. 
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2.1 


THE  e  DETECTION  STATISTIC 

If  a  simple  energy  detector  were  used  with  the  array 
output  to  detect  a  possible  test,  a  time-varying  threshold  would 
be  necessary  to  take  advantage  of  the  general  decay  in  the 
energy  level  of  the  interfering  earthquake  signal.  As  a  simpler 
approach,  we  have  considered  an  energy  ratio  statistic  and  a 
constant  threshold.  This  detection  statistic  is  defined  as  the 
ratio  of  the  average  power  in  a  one-second  interval  to  that  of 
the  previous  five-second  interval.  The  five-second  integration 
interval  was  chosen  as  adequate  to  smooth  out  the  short  term 
fluctuations  that  are  apparent  in  the  energy  decay  curves. 

The  e  statistic  is  calculated,  as  a  function  of  time, 
from  the  array  output.  The  array  output  is  calculated  for  each 
of  several  possible  test  epicenters  by  aligning  the  seismograms 
so  that  signals  from  the  test  epicenter  would  arrive  in  phase, 
normalizing  the  seismograms,  and  adding  them.  Two  normalizations 
of  the  seismograms  have  been  used.  For  the  first,  all  seismograms 
are  normalized  to  have  the  same  P-phase  energy  (defined  by  a  one- 
second  time  window).  In  the  case  where  some  stations  have  poor 
signal-to-noise  ratios,  this  normalization  can  lead  to  excessive 
noise  in  the  array  output.  For  times  well  after  the  earthquake, 
when  the  coda  energy  is  of  the  same  order  as  the  noise,  this 
noise  is  a  significant  limitation  on  detection  performance.  This 
problem  can  be  alleviated  by  using  maximal-ratio  combining,  in 
which  each  record  is  weighted  by  the  ratio  of  the  square  root  of 
the  P-phase  energy  to  the  average  power  in  the  noise  preceding  the 
P-phase.  One  potential  problem  in  using  the  maximal-ratio  combining 
is  that  the  directional  properties  of  the  array  may  be  degraded. 

As  will  be  illustrated  in  the  data  presented  below,  it  appears 
that  the  maximal-ratio  combining  is  to  be  preferred,  despite  this 
potential  disadvantage. 

Figure  2.1  illustrates  the  behavior  of  e  as  a  function  of 
time  when  the  array  is  steered  at  the  earthquake's  epicenter;  this 
figure  is  based  on  maximal-ratio  combining.  The  maximum  excursion 
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is  1700,  and  the  only  other  large  excursion  occurs  roughly 
15  seconds  after  the  arrival  of  the  event. 

The  same  curve  has  been  calculated  with  the  array 
steered  to  each  of  thirty  grid  points  surrounding  the  earth¬ 
quake's  epicenter.  These  grid  points  were  defined  by  six 
bearings  and  distances  of  1,  3>  5}  7  and  9  degrees  from  the 
earthquake's  epicenter.  The  results  of  these  calculations 
indicate  average  levels  of  e  very  close  to  1.0,  except  in 
time  intervals  corresponding  to  arrivals  of  some  phases  from 
the  earthquake.  The  rms  value  of  this  statistic  is  also  close 
to  unity. 

In  order  to  set  a  threshold  that  excludes  all  but  a 
few  false  alarms  it  is  necessary  to  tabulate  the  maximum 
excursions  of  e  as  a  function  of  the  test  epicenter  to  which 
the  array  is  steered.  Figure  2.2k  gives  these  maximum  values 
at  each  grid  point  surrounding  the  earthquake's  epicenter  for 
the  case  of  equal  P-phase  energy  normalization.  Figure  2.2B 
is  the  same  presentation  based  on  maximal-ratio  combining. 
Maximum  excursions  as  high  as  50  in  Figure  2.2k  and  150  in 
Figure  2.2B  occur. 

These  maximum  excursions,  and  hence  the  allowable 
detection  thresholds,  may  be  reduced  significantly  by  excluding 
the  times  corresponding  to  P-phase  arrivals  from  these  calcu¬ 
lations.  As  discussed  in  more  detail  in  a  previous  report  [3] 
this  region  of  time  can  be  more  effectively  monitored  by  a 
network  of  single  stations  surrounding  the  region  of  interest. 
For  this  reason,  we  shall  exclude  these  times  from  our  calcu¬ 
lations  of  array  characteristics.  Excluding  those  times  that 
could  be  monitored  by  a  single  station  network  that  is  well 
distributed  in  azimuth  and  has  an  average  range  of  ^0°  from  the 
earthquake's  epicenter,  leads  to  the  maximum  e  values  given 
in  Figure  2.3. 

In  Figures  2.2  and  2.3,  the  maximal-ratio  combining 
leads  to  larger  maximum  excursions  of  e  than  the  equal  P-phase 
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FIGURE  2.2  B 

MAXIMUM  EXCURSIONS  OF  €  STATISTIC 
INCLUDING  P-PHASE  ARRIVALS 


-10- 


FIGURE  2.3  A 

MAXIMUM  EXCURSIONS  OF  €  STATISTIC 
EXCLUDING  P-PHASE  ARRIVALS 


-11- 


7 


FIGURE  2.3  B 
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energy  normalizations.  If  we  wanted  to  allow  a  margin  of  safety 
of  a  factor  of  two,  we  would  have  to  choose  a  detection  threshold 
of  at  least  26  for  the  equal  P-phase  energy  case  and  36  for  the 
maximal-ratio  case.  Since  these  numbers  seem  rather  large,  we 
have  also  considered  an  alternate  detection  statistic,  n ,  in  the 
hope  that  it  would  turn  out  to  be  a  more  sensitive  statistic. 

2.2  THE  n  DETECTION  STATISTIC 

A  second  detection  statistic,  n ,  is  defined  as  the 
product  of  c  and  the  average  paired  correlation  coefficient 
calculated  from  the  corresponding  one-second  interval  of  all  of 
the  seismograms.  The  objective  in  using  this  statistic  is  the 
reduction  of  "spurious"  peaks  in  the  array  output  compared  to 
"legitimate"  peaks  that  result  from  correlated  arrivals  at  each 
station . 

Some  properties  of  the  n  statistic  may  be  observed  from 
Figure  2.4,  which  is  a  plot  of  the  logarithm  of  the  magnitude 
of  n  vs.  time.  Since  maximal-ratio  combining  is  used  in  both 
cases,  each  value  of  n  in  Figure  2.4  is  just  the  product  of  the 
corresponding  value  of  e  in  Figure  2.1  and  the  average  paired 
correlation  coefficient  for  that  time.  Figure  2.4  indicates 
that  the  n  statistic  is  generally  quite  low  compared  to  the  e 
statistic  except  at  the  correlated  P-phase  arrival  (recall  that 
the  beam  is  steered  to  the  earthquake  epicenter)  and  at  the 
secondary  peak  15  seconds  later.  (Based  on  earlier  studies[4]'  of 
this  event,  we  do  not  think  that  this  secondary  peak  is  pP 
since  pP  has  been  found  within  five  seconds  of  P.) 

Figure  2.5  shows  the  maximum  of  the  n  detection  sta¬ 
tistic  when  the  total  records  are  considered  (including  the  P- 
phase  arrivals).  It  is  apparent  that  these  values  are  signi¬ 
ficantly  smaller  than  the  e's  shown  in  Figure  2.2.  There  is 
still  a  large  excursion  of  n  in  Figure  2 . 5B  at  the  grid  point 
Az=0  and  R=7 .  At  this  same  grid  point  the  corresponding  maximum 
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value  of  e  was  found  to  be  133.  Figure  2.6  shows  the  maximum 
excursions  of  this  statistic  when  the  time  interval  that  would  be 
covered  by  the  single-station  network  is  excluded.  The  maximum 
excursions  obtained  indicate  that  in  the  case  of  maximal-ratio 
combining,  an  n  of  2.0  (3  dB  greater  than  the  maximum)  is  a 
conservative  detection  threshold.  Using  the  equal  P-phase 
energy  normalization,  the  corresponding  threshold  would  be  6.0. 

2.3  COMPARISON  OF  THE  PERFORMANCE  OF  THE  e  AND  n 

DETECTION  STATISTICS 

The  preceding  figures  give  some  indication  of  what 
thresholds  could  be  used  with  e  and  with  n  in  order  to  avoid 
false  alarms.  This  is,  of  course,  only  part  of  the  story.  We 
now  wish  to  examine  the  sensitivity  of  these  two  statistics  to 
the  actual  presence  of  a  second  seismic  event  in  the  coda  of  the 
first.  To  study  this  question  several  simulations  have  been  per¬ 
formed.  In  these  simulations  a  scaled  replica  of  the  P-phase  of 
each  seismogram  is  added  to  the  same  seismogram  at  a  time  corres¬ 
ponding  to  the  arrival  of  a  hypothetical  test  that  is  set  off 
60  seconds  after  the  earthquake  and  is  located  at  a  designated 
grid  point . 

The  results  of  these  calculations  are  summarized  in 
Figures  2.7  and  2.8.  Each  graph  contains  six  curves,  one  for 
each  grid  point  at  a  fixed  range.  Each  curve  represents  the  locus 
of  pairs  of  values  of  e  and  n  as  the  scale  factor  for  the  added 
P-phase  is  varied.  Thus  these  plots  allow  us  to  examine  the  rela¬ 
tive  values  of  e  and  n  as  the  magnitude  of  the  (simulated)  second 
event  is  increased.  It  may  be  observed  from  these  figures  that 
the  general  form  of  the  e-n  relation  does  not  vary  a  great  deal 
from  grid  point  to  grid  point.  Also  indicated  on  these  plots 
are  the  thresholds  suggested  by  the  earlier  figures  as  sufficient 
to  avoid  false  alarms.  It  is  clear  from  these  plots  that  the  n 
statistic  will  cross  its  "false-alarm  proof"  threshold  before  the 
e  statistic  crosses  its  correspondingly  conservative  threshold. 
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FIGURE  2.6  B 
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FIGURE  2 .7  B 
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(RECORDS  NORMALIZED  TO  SAME  P-PHASE  ENERGY) 


-21- 


€ 


FIGURE  2.7 C 

RELATION  BETWEEN  €  AND  T) 

(RECORDS  NORMALIZED  TO  SAME  P- PHASE  ENERGY) 


-22- 


€ 


FIGURE  2.8 A 
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The  difference  is  only  minor  in  the  case  of  the  equal  P-phase 
energy  normalization,  but  it  is  quite  large  in  the  case  of  the 
maximal-ratio  combining.  For  example,  in  Figures  2 . 8A  and  B, 
at  the  point  where  the  n  threshold  is  exceeded,  e  is  at  only 
about  one-fifth  of  its  threshold.  Therefore,  comparing  both 
statistics  on  the  basis  of  a  conservative  (no  false  alarms) 
threshold,  it  appears  that  the  n  statistic  is  more  sensitive 
by  about  a  factor  of  five  in  power  (roughly  7  dB).  Furthermore, 
at  the  level  at  which  detection  occurs,  e  is  less  than  10. 

This  implies  that  detection  is  assured  if  a  one-second  interval 
of  the  array  output  has  a  power  level  that  is  10  dB  above  that 
of  the  preceding  five  seconds,  provided  that  the  maximal-ratio 
array  output  is  used,  and  the  n  statistic  is  used  for  detection. 
With  this  observation  we  can  approximately  describe  the  sensi¬ 
tivity  of  the  array,  as  a  function  of  grid  point  and  delay, 
by  examining  the  average  energy  in  a  five-second  time  window 
as  a  function  of  time  after  the  earthquake  for  each  grid  point. 

For  a  given  time  and  grid  point  it  is  assumed  that  the  minimum 
detectable  shot  energy  is  10  dB  larger  than  the  corresponding 
five-second  average  energy, 

2  A  MINIMUM  DETECTABLE  MAGNITUDE  AS  A  FUNCTION  OF  DELAY 

TIME  AND  GRID  POINT 

As  suggested  by  the  preceding  paragraph  the  array's 
detection  sensitivity  for  each  grid  point  as  a  function  of 
time  can  be  calculated  by  computing  average  power  curves  for 
the  array  output  that  are  smoothed  over  five  seconds,  and  adding 
the  10  dB  required  for  detection  to  these  curves  Some  results 
of  this  calculation  are  presented  in  Figures  2.9  through  2.11. 
Figure  2.9  shows  the  detectable  magnitude  of  a  nuclear  test  when 
the  test  site  coincides  with  the  earthquake's  epicenter.  From 
this  figure  it  can  be  observed  that  a  nuclear  explosion  of  the 
same  magnitude  as  the  earthquake  would  not  be  detectable  if  the 
test  site  and  the  earthquake  epicenter  coincided  and  the  explosion 
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DETECTABLE  MAGNITUDE  VS.  DELAY 
(TEST-SITE  LOCATED  AT  EPICENTER) 
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DETECTABLE  MAGNITUDE  VS.  DELAY 
(TEST-SITE  LOCATED  AT  AZ  =  240°,  R  =  3°) 
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(TEST-SITE  LOCATED  AT  AZ  =  0°,  R=5°) 


occured  within  5-10  seconds  after  the  earthquake.  Similarly, 
it  can  be  observed  that  a  nuclear  test  which  is  0.5  magnitude 
units  smaller  than  the  natural  event  will  be  detected  if  it 
occurs  later  than  50  seconds  after  the  natural  event  and  the 
epicenters  coincide. 

This  same  calculation  performed  for  all  of  the  grid 
points  considered  above  indicates  that  a  nuclear  test  with  the 
same  magnitude  as  the  natural  event  considered,  could  be  detected 
at  all  of  the  other  grid  points,  independent  of  the  earthquake- 
shot  delay.  Figures  2.10  and  2.11  show  the  detectable  magnitude 
as  a  function  of  delay  for  two  representative  grid-points, 

Az  =  0° ,  R  =  5°  and  Az  =  2^0°,  R  =  3° . 

2.5  TIME  INTERVALS  AVAILABLE  FOR  CONCEALING  NUCLEAR  TESTS 

The  monitoring  capabilities  of  the  continental-size 
array  will  now  be  summarized  in  terms  of  the  time  intervals  over 
which  a  nuclear  test  can  be  detonated  without  risk  of  detection 
for  various  points  in  the  epicenter  region.  Figure  2.12  shows  the 
time-intervals  outside  of  which  the  detection  system  formed  by 
a  continental-size  array  and  a  single  station  network  will  detect 
an  explosion  of  magnitude  5.0  for  different  test  site  locations. 

As  in  earlier  figures,  the  earthquake's  epicenter  is  in  the 
center  of  the  figure.  The  two  numbers  shown  at  each  grid  point 
correspond  to  the  minimum  delay  imposed  by  the  single-station 
network  and  the  maximum  delay  imposed  by  the  continental-size 
array.  The  minimum  waiting  time  indicated  in  Figure  2.12  is  a 
consequence  of  the  assumed  monitoring  by  a  single  station  net¬ 
work,  that  is,  well  distributed  in  azimuth  and  located  at  an 
average  range  of  ^0°  from  the  epicenter.  Grid  points  marked 
with  an  X  identify  the  test  site  locations  for  which  the  minimum 
waiting  time  is  greater  than  the  maximum  imposed  by  the  array. 
Observe  that  the  system  considered  will  detect  a  5.0  magnitude 
explosion  occurring  with  a  delay  greater  than  50  seconds.  When 
the  test  site  does  not  coincide  with  the  epicenter  of  the  natural 
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TIME  INTERVALS  AVAILABLE  FOR  CONCEALING 
(  TEST  MAGNITUDE  5.0  ) 
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event,  the  time  intervals  in  which  explosions  of  5.0  magnitude 
would  not  be  detected  are  shorter  than  50  seconds . 

The  continental-size  array  that  is  assumed  in  these 
calculations  is  located  to  the  east  of  the  earthquake's  epi¬ 
center.  Therefore,  grid  points  shown  on  the  right-hand  side 
of  Figure  2.12  are  closer  to  the  array  than  the  earthquake 
epicenter.  It  can  be  observed  that  it  is  easier  to  detect  an 
explosion  when  the  array  and  the  test  site  are  located  on  oppo¬ 
site  sides  of  the  epicenter  of  the  earthquake,  than  otherwise. 

The  presentation  of  Figure  2.12  is  repeated  in  Figure 
2.13  for  the  case  of  a  mangitude  4.5  explosion,  i.e.,  one  magni¬ 
tude  unit  less  than  the  natural  event.  The  minimum  delay 
imposed  by  the  single  station  network  is,  of  course,  the  same 
as  in  Figure  2.12.  However,  the  maximum  allowable  delay  times 
are  considerably  larger.  It  can  be  seen  that,  for  coincident 
epicenter  and  test  site,  a  nuclear  explosion  would  be  detected 
with  the  system  described  only  if  it  occurred  with  a  delay 
greater  than  about  three  minutes. 

An  additional  constraint  on  the  maximum  earthquake- 
nuclear  explosion  delay  could  be  derived  by  considering  the 
single-seismometer  stations  that  are  located  on  the  opposite 
side  of  the  epicenter  from  the  continental  array.  Signals  arriv¬ 
ing  at  these  stations  from  grid  points  situated  between  the 
epicenter  and  the  array  might  arrive  with  a  considerable  delay. 

As  an  example,  a  nuclear  explosion  of  magnitude  4.5  and  which 
occurs  at  a  site  located  on  the  grid  point  Az  =  60°  and  R  =  9° 
with  a  delay  of  260  seconds,  will  arrive  at  some  of  the  single 
stations  (those  for  which  the  azimuth  from  the  epicenter  is 
approximately  240°)  with  a  delay  of  approximately  330  seconds. 
Delays  of  this  size  might  lead  to  detection  of  the  signal  arriving 
from  the  nuclear  test  on  the  record  obtained  from  the  single 
seismometer . 
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2.6  SUMMARY  AND  CONCLUSIONS 

The  capabilities  of  a  hypothetical  detection  network, 
consisting  of  a  continental-size  array  and  a  network  of  single 
stations,  have  been  summarized  in  terms  of  the  allowable  time 
intervals  during  which  a  nuclear  test  could  be  detonated  in  the 
vicinity  of  an  earthquake  without  risking  detection,  as  a  function 
of  the  magnitude  of  the  test  and  the  position  of  the  test  site 
relative  to  the  array  and  the  earthquake  epicenter.  These  inter¬ 
vals  are  presented  in  Figures  2.12  and  2.13,  for  shots  that  are 
1/2  and  1  magnitude  unit  smaller  than  the  earthquake,  The  limita¬ 
tions  imposed  by  these  intervals  suggest  that  the  practical  prob¬ 
lems  of  concealing  a  nuclear  test  from  a  detection  network  of  the 
sort  described  here  would  be  quite  formidable. 

Unfortunately,  an  enormous  number  of  computations  would 
be  necessary  to  implement  the  detection  network.  If  it  were  only 
necessary  to  monitor  a  small  number  of  potential  test  sites,  the 
computation  problem  would,  of  course,  be  more  manageable.  One 
possibility  for  reducing  the  number  of  computations  would  be  to 
increase  the  beamwidth  of  the  array  by  reducing  its  aperture,  but 
this  would  presumably  degrade  the  sensitivity  of  the  system. 

The  calculations  presented  in  this  section  are  all  based 
on  data  from  a  single  earthquake.  As  far  as  we  know,  there  is 
nothing  unusual  about  this  earthquake  and  these  results  are 
representative  of  what  could  be  obtained  with  other  earthquakes. 

It  would  probably  be  useful,  however,  to  process  data  from  another 
event  in  order  to  confirm  this  conjecture  and  to  see  how  much 
detailed  variation  there  is  between  events.  We  do  not  anticipate 
that  the  detailed  fluctuations  between  events  would,  in  themselves, 
be  of  interest,  except  that  knowledge  of  their  extent  would  allow 
conservative  detection  levels  to  be  estimated  with  some  confidence 
that  they  would  be  generally  applicable. 
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SECTION  III 


BEAMFORMING  FOR  EPICENTER  LOCATION  -  THEORETICAL  BACKGROUND 

The  problem  of  estimating  epicenters  of  seismic  events 
using  only  data  from  LASA-Montana  has  been  studied  by  several 
groups.  In  an  earlier  report  [5  ] ,  we  presented  experimental 
results  based  on  a  beamsplitting  technique  and  some  tentative 
conclusions  concerning  the  dependence  of  the  precision  of  the 
epicenter  estimates  on  certain  parameters.  As  discussed  in  that 
report,  the  original  motive  for  the  beamsplitting  technique  was 
that  the  coherent  nature  of  the  processing  would  yield  better 
results  (perhaps  at  the  expense  of  additional  complexity)  than 
could  be  obtained  with  conventional  plane  wave  fits  that  were 
based  upon  independent  time  picks  on  the  separate  (noisy)  seismo¬ 
grams.  In  many  cases,  at  least,  it  now  appears  that  travel-time 
anomalies,  not  additive  noise,  are  the  major  contributor  to  errors 
in  epicenter  location.  In  these  cases  the  original  motive  for 
using  the  beamsplit  technique  is  no  longer  obviously  valid.  With 
this  in  mind,  as  well  as  our  earlier  experience  reported  previ¬ 
ously,  we  have  been  re-examining  the  beamsplitting  technique  for 
locating  epicenters  with  a  view  toward  learning  the  irreducible 
limits  on  location  precision  and  simplifying  the  complexity  of  the 
calculations  used  to  make  these  estimates.  In  this  section,  we 
shall  present  theoretical  analyses  of  some  rather  idealized 
problems  that  are  related  to  the  epicenter  location  problem. 

Results  of  simulation  studies  are  discussed  in  Section  IV. 

The  contents  of  this  section  are  organized  as  follows. 
Section  3.1  states  the  assumptions  that  are  common  to  all  of  the 
derivations.  The  most  important  of  these  are  that  an  identical 
signal  waveform  is  present  at  each  seismometer,  the  noise  waveforms 
are  independent  but  share  a  common  statistical  description,  and 
the  signal-to-noise  ratio  is  large  enough  that  all  errors  are 
small.  Section  3-2  summarizes  the  form  of  a  least  squares 
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(regression)  fit  to  a  set  of  time  picks,  and  gives  the  variance 
in  the  estimate  of  a  parameter,  m,  that  is  proportional  to  the 
cosine  of  the  angle  of  incidence,  in  terms  of  the  variance  of  the 
estimates  of  the  individual  time  picks.  Section  3-3  summarizes 
the  performance  of  three  methods  of  estimating  individual  arrival 
times:  a  zero  crossing  pick,  a  maximum  value  pick,  and  a  cross¬ 

correlation  with  a  perfect  replica  of  the  signal.  Section  3D 
discusses  the  beamforming  calculation  and  demonstrates  its  equi¬ 
valence  to  the  plane  wave  fit  following  a  correlation  time  pick. 

The  important  point  here  is  that  the  beamforming  technique  does 
not  require  prior  knowledge  of  the  detailed  waveshape  of  the  signal, 
whereas  the  correlation  time  pick  does  require  this  knowledge. 
Finally,  Section  3-5  summarizes  the  results  in  terms  of  a  few 
examples.  In  one  of  these,  a  21-element,  200-km  array  with  a 
20  km/sec  phase  velocity,  and  a  10  dB  signal-to-noise  ratio, 
leads  to  a  standard  deviation  in  angular  measurement  of  0.3°. 

The  details  of  the  more  involved  derivations  presented  in  this 
section  are  presented  in  Appendices  A  and  B. 

3.1  STATEMENT  OF  ASSUMPTIONS 

We  assume  a  uniform  linear  array  with  N  elements  indexed 
by  i.  The  waveform  at  the  i^  element  is  assumed  to  be  given  by 

x1(t)  =  s(t-im+T^)+n^( t)  (3-D 

where  the  n^(t)  are  independent  stationary  noise  processes 
with  a  common  power  density  spectrum,  S  (f),  the  t.  are  inde¬ 
pendent  ,  identically  distributed  random  variables  that  correspond 
to  perturbations  on  the  arriving  wavefront,  and  s(t)  is  given  by 

s ( t )  =  sin  2nf  t  -T/2  <  t  <  T/2  (3-2) 

=  0  elsewhere 

where  T  is  the  period  of  the  sinewave.  The  parameter  m  in 
Equation  (3-1)  is  related  to  the  velocity,  c,  the 
element  spacing,  A,  and  the  angle  of  incidence,  6,  by 

sine  =  m  (3-3) 

where  6  and  c  are  both  measured  in  the  plane  of  the  array. 
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For  the  application  of  interest  here,  fg  -  1  Hz,  but  the  parameter 
will  be  retained  in  the  formulation  in  order  to  show  parameter 
dependence  and  allow  scaling  of  results. 

We  shall  consider  two  general  schemes  for  estimating 
the  parameter  m  in  Equation  (3-1).  The  first  will  consist  of 
estimating  the  "arrival  time,"  y^,  from  each  of  the  x^(t);  in 
the  absence  of  noise  this  estimate  would  be 

y±  =  im--^  (3-4) 

The  are  then  combined,  e.g.,  by  a  least  squares  regression, 
to  estimate  m.  (It  is  for  this  reason  that  we  shall  find  it 
convenient  to  refer  to  m  as  the  slope.  It  is  the  slope  of  the 
correct  regression  line;  it  is  not  the  slope  of  the  wavefront.) 
The  second  general  scheme  consists  of  forming  a  delayed  sum  for 
the  array ,  g( t ,m) . 

N 

g(t,m)  =  l  x  (t+mi)  (3-5) 

i  =  1  1 

~  p 

and  then  choosing  m  to  maximize  the  integral  of  g  (t,m)  over 

some  time  interval - 


3.2  SLOPE  ESTIMATE  FROM  UNRELIABLE  TIME  PICKS 

It  is  assumed  that  the  observed  arrival  times  are  given  by 


y_.  =  im+b  +  z^ 


(3-6) 


where  the  z^  are  zero-mean  random  variables  with 


z  .  z  .  =  6  .  , 
i  J  ij 

and  m  and  b  are  unknown  constants. 

One  way  to  estimate  m  is  to  choose  in  and  b  in  such 
way  as  to  minimize  the  mean  square  error,  e(m,b): 


(3-7) 


a 


£ (m,b) 


1 

N 


N 


l  (y.-im-b)' 
i=l  1 


(3-8) 


This  defines  the  standard  regression  fit  to  the  data  points, and 
it  is  easy  to  show  that  m  is  given  by  the  following  expression. 
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m  = 


?  n  N+1^ 

1  (i  — p-)yi 

i  =  l _ d 

Yj  n(n2-i) 


(3-9) 


which  may  be  rewritten  using  Equation  (3-6)  to  yield 

N  N+l  ■ 


m  =  m  + 


I  a  - 

i  =  l _  1 

j4  N(N  -1) 


(3-10) 


The  variance  of  m  may  be  calculated  from  Equation  (3-10)  and 
the  result  is 

12 


var[m]  = 


N(N  -1) 


(3-11) 


It  should  be  recalled  that  this  is  based  on  the  z^  having  unity 
variance;  Equation  (3-11)  may  be  generalized  by  multiplying 
the  right  side  by  the  variance  of  the  z's. 

Several  observations  about  the  estimate  given  by 
Equation  (3-10)  may  be  made.  The  estimate  was  derived  by 
minimizing  the  mean  square  error  given  by  Equation  (3-8).  It 
turns  out  that  the  estimate  weights  the  time  measurements 
linearly  with  distance  from  the  center  of  the  array,  as  would 
seem  intuitively  reasonable.  In  the  case  that  the  z^  are 
Gaussian, the  estimate  given  by  Equation  (3-9)  may  be  shown  to 
be  maximum  likelihood  and  minimum  variance.  In  the  case  where 
the  z^  are  not  Gaussian,  it  is  not  obvious  that  the  least 
squares  criterion  results  in  a  minimum  variance  estimate  of  m. 

One  generalized  class  of  estimates  can  be  obtained  by  changing 
the  weights  applied  to  the  y^  in  Equation  (3-9)  (and  appro¬ 
priately  correcting  the  denominator  to  maintain  an  unbiased 
estimator) .  It  turns  out  that  this  is  equivalent  to  mini¬ 
mizing  a  weighted  error  function: 

1  N  P 

ew(™>6)  =  jj-  l  (y1-im-b)^ai  (3-12) 
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It  may  be  shown  that  of  all  possible  choices  of  the  weights,  a^ , 
the  uniform  weights  (a.  =  1)  lead  to  the  minimum  variance  esti¬ 
mate  of  m  [6].  Thus  the  estimate  of  Equation  (3-9)  appears 
reasonable  for  more  than  just  Gauss  errors. 


3.3  ACCURACY  OF  TIME  PICKS  ON  INDIVIDUAL  RECORDS 

There  are  several  ways  In  which  the  arrival  time  for 

a  given  record  might  be  estimated.  In  this  section  we  shall 

discuss  the  statistical  properties  of  three  methods.  Our 

objective  in  considering  these  methods  is  to  indicate  how 

the  variance  of  the  time  picks  (and  hence  that  of  m)  might 

depend  on  the  level  and  shape  of  the  noise  spectrum.  As  before 

(see  Equations  3-1  and  3-2)  we  assume  that  each  record  consists 

of  one  period  of  a  sine  wave  plus  stationary  noise, and  it  is 

desired  to  estimate  the  "arrival  time"  of  the  signal.  The 

three  cases  derived  in  Appendix  A  are:l)  the  zero  crossing 

near  the  center  of  the  transient  is  selected,  2)  the  point  of 

maximum  value  near  the  end  of  the  transient  is  selected,  and 

3)  the  record  is  correlated  with  a  perfect  replica  of  s(t)  and 

the  time  of  maximum  correlation  is  selected.  The  variances  of 

2  2  2 

the  time  picks  due  to  the  noise  are  given  as  o^,  o 2,  and  o.^ 
for  the  three  methods.  (See  Equations  A-8,  A-9  and  A-17  of 
Appendix  A.)  In  all  three  derivations,  large  signal-to-noise 
ratios  and  small  errors  are  assumed. 

o\  =  jSn(f)df  (3-13) 

“0 


2 

0  2 


2 

°3 


4  Aft>2  vf)df 

Wq  0 


/[|(sinir 

u0 


J-)  (  n°-,)]2s  (f)df 

0  fg-f 


(3-1^) 

(3-15) 


where  u  =  2 irf  is  used  for  convenience  in  some  of  the  expressions. 
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A  quick  check  on  the  consistency  of  Equations  ( 3-13 )-( 3-15 ) 
is  the  case  of  extremely  narrowband  noise,  for  which  Sn(f)  consists 
only  of  impulses  at  if,.  In  this  case  all  three  variances  are 
equal,  as  would  be  expected.  It  does  appear  to  be  possible  to 
choose  a  power  density  spectrum  for  the  noise  such  that  the  third 
case  is  worse  than  the  first  case,  but  for  most  reasonable  assump¬ 
tions  about  the  power  density  spectrum  of  the  noise,  one  would 
guess  that  the  correlation  technique  would  lead  to  a  better  esti¬ 
mate  than  the  zero  crossing  estimate.  This  seems  like  a  reasonable 
intuitive  guess  simply  because  the  correlation  estimate  involves 
an  integration  over  time  that  should  help  to  reduce  the  effects  of 
the  noise. 

It  is  difficult  to  make  general  conclusions  about  the 
relative  size  of  the  three  variances  without  making  specific 
assumptions  about  S  (f),  since  the  multipliers  of  Sn(f)  in  each 

of  the  integrands  of  Equations  ( 3— 1 3 ) — ( 3-15)  do  not  satisfy  the 

2  2 

same  rank  ordering  for  all  f.  Comparing  and  it  is  clear 

that  the  shape  of  S  (f)  is  important.  The  time  pick  based  on  the 
maximum  is  less  (more)  sensitive  to  noise  below  (above)  f  q  than 
the  time  pick  based  on  the  zero  crossing.  The  reasonableness  of 
this  result  can  be  understood  by  considering  the  effects  of  the 

extremes  of  a  dc  offset  or  a  low  level  high  frequency  noise  on  the 

2  2 

two  picks.  Comparing  and  is  a  little  more  difficult.  The 
multiplier  in  the  integrand  of  Equation  (3-15)  is  less  than  unity 
for  most,  but  not  all,  f.  (It  equals  unity  at  f=fQ  and  reaches 

a  maximum  at  a  slightly  higher  frequency.)  Therefore,  it  would 

2  2 

be  expected  in  general  that  but  there  do  exist  (rather 

artificial)  cases  where  the  reverse  is  true. 

If  in  the  third  case,  white  noise  of  spectral  level  Nq  , 
is  assumed, 

Sn(f)  =  Nq  (3-16) 

The  variance  of  the  time  estimate  turns  out  to  be 

=  N0/™0  (3-17) 
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In  order  to  compare  this  result  to  the  variance  that  would 
obtain  in  the  case  of  a  zero  crossing  time  pick,  assume  that 
the  power  density  spectrum  of  the  noise  is  given  by 


Sn(f)  =  NQ  | f |  <  B 

=  0  |  f  |  >  B 

Then  it  follows  that 

N 

2  0  2  ttB 

O  =  -  - 


(3-18) 


(3-19) 


The  motive  for  the  different  assumptions  about  the  power 
density  spectrum  is  that  the  correlation  technique  implies  a 
filtering  of  its  own,  whereas  the  zero  crossing  time  pick 
does  not.  Clearly,  if  there  were  a  good  deal  of  noise  energy 
in  a  different  region  of  the  spectrum  than  the  signal  energy, 
it  would  be  sensible  to  fiiter  it  out  before  choosing  the 
time  of  the  zero  crossing.  The  parameter  B  in  Equation  (3-18) 
would  presumably  be  somewhat  larger  than  f q ,  but  not  very 
much  larger.  This  being  the  case,  it  may  be  seen  that  the 
variance  for  the  zero  crossing  pick  is  not  very  much  larger 
than  that  of  the  correlation  pick.'  Furthermore,  in  writing 
Equation  (3-18)  we  have  excluded  the  possibility  of  filtering 
out  the  very  low  frequency  region  of  the  noise,  which  presumably 
could  be  done  without  serious  distortion  of  the  signal.  The 
principal  conclusion  intended  from  Equations  (3-17)  and  (3-19) 
is  that  while  the  correlation  time  pick  is  usually  superior 
to  that  of  the  zero  crossing  time  pick,  it  appears  that  the 
principal  reason  for  this  is  that  a  narrowband  filtering  is 
involved,  and  if  an  equivalent  narrowband  filtering  is  inserted 
before  the  zero  crossing  time  pick,  relatively  little  difference 
would  exist  between  these  two  estimates.  It  should  be  stressed 
that  all  of  these  comments  are  based  on  the  assumption  of 
identical  signals  at  each  element.  If  this  is  not  the  case,  the 
correlation  pick  would  presumably  be  degraded.  Furthermore,  in 
the  case  of  the  zero  crossing  pick  with  nonidentical  signals,  the 
effect  of  prefiltering  would  have  to  be  carefully  considered. 


3 . 4  BEAMFORMING 

The  beamforming  technique  of  estimating  epicenters  for 
seismic  events  is  fundamentally  different  from  the  plane-wave- 
fit  technique  discussed  above.  In  the  beamforming  technique  the 
array  is  "steered"  to  several  trial  epicenters  and  the  epicenter 
yielding  the  largest  energy  in  the  array  output  is  chosen  as  the 
estimate  of  the  true  epicenter.  The  steering  of  the  array  is 
accomplished  by  aligning  each  seismogram  before  addition  in  such 
a  way  that  if  the  test  epicenter  were  the  true  epicenter,  the 
signal  portion  of  each  waveform  would  add  coherently  in  the  array 
output.  In  general,  it  would  be  expected  that  the  beamforming 
technique  involves  more  computations  than  the  conceptually  simpler 
plane-wave  fit  based  on  separate  time  picks  from  each  channel. 

More  calculations  are  involved  because,  unless  some  simplification 
is  developed,  the  array  output  must  be  calculated  for  several 
steerings  of  the  beam  and  each  array  output  is  based  on  a  delay- 
and-sum  computation.  The  question  naturally  arises  as  to  why 
beamforming  should  be  used  for  the  location  of  seismic  events. 

The  principal  reason  is  that  it  may  be  a  useful  technique  for 
reducing  the  noise  contribution  to  epicenter  location  errors. 

This  would  be  expected  because  the  beamforming  technique  is 
coherent  in  the  sense  that  all  of  the  signals  are  added  prior  to 
any  estimation  procedures  and,  in  general,  predecision  processing 
is  preferable  to  other  alternatives.  A  secondary  reason  for 
using  beamforming  rather  than  the  plane-wave  fit  is  that  the  best 
of  the  time  pick  procedures  discussed  above  depends  upon  perfect 
knowledge  of  the  waveshape  of  the  signal  component  in  each  record. 
The  beamforming  technique  does  not  require  perfect  knowledge  of 
the  signal,  although,  in  its  simplest  form,  it  does  assume  that 
the  same  signal  is  present  at  each  seismometer.  It  will  be  shown 
below  that  in  the  case  of  small  errors  the  beamforming  technique 
is  essentially  equivalent  in  performance  to  the  correlation  time 
pick  followed  by  a  plane-wave  fit.  This  derivation  will  be 
presented  to  indicate  the  potential  value  of  the  beamforming 
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technique:  results  comparable  to  those  obtainable  from  the  best 

of  the  other  techniques  may  be  achieved  without  prior  knowledge 
of  the  signal  waveform 

In  the  two-element  case  the  beamforming  estimate  of  the 
slope,  m,  is  obtained  by  maximizing  the  following  expression: 

/ C x ^ ( t )  +  x2 ( t+m ) 1 “~dt  =  J[s(t)+n  ( t )  +  s ( t  +  m-m)  +  np ( t  +  m) ]^dt  (3-20) 

In  this  derivation  it  is  assumed  that  the  limits  of  the  integral 
loosely  bound  the  duration  of  the  signal  term.  In  this  deriva¬ 
tion,  as  well  as  all  others  in  this  section,  a  large  signal-to- 
noise  ratio  will  be  assumed  so  that  there  is  little  difficulty 
in  approximately  locating  the  signal  and  thereby  setting  the 
appropriate  limits  on  the  integral  The  choice  of  m  in  Equa¬ 
tion  (3-20)  influences  the  value  of  the  integral  only  through 
the  cross-terms  obtained  by  expanding  the  square.  Therefore, 
Equation  (3-20)  will  be  maximized  when  the  expression  below  is 
maximized . 

/[  s  ( t )  s  ( t+m-m)  +  s  ( t )  n2  ( t+m  )+s  ( t+m  )n-^(  t  )+n ,( t  )n2  ( t+m)  ]dt  (3-21) 

If  we  neglect  the  second  order  noise  term,  what  remains  is  an 
expression  very  similar  to  the  expression  that  would  appear  in 
the  case  of  a  time  pick  via  correlation  with  a  perfect  replica 
of  the  signal.  In  the  case  of  the  correlation  time  pick  (see 
Appendix  A)  there  is  only  a  single  signal-noise  term  Instead  of 
the  two  cross-terms  that  appear  in  expression  (3-21).  Otherwise 
the  expression  given  in  (3-21)  is  identical  to  the  one  that  appears 
in  the  case  of  a  correlation  time  pick  From  this  observation  it 
follows  that  the  variance  of  the  estimate  m  in  this  problem  will 
be  twice  the  variance  calculated  for  the  correlation  time  pick 
and  given  by  Equation  (3-15).  If  the  slope  m  were  estimated  by 
using  the  correlation  time  pick  on  each  of  the  two  elements  and 
then  drawing  a  straight  line  between  the  resulting  estimates,  the 
variance  in  the  resulting  estimate  would  also  be  twice  that  of 
the  estimate  of  the  individual  time  picks.  Therefore,  it  may  be 
concluded  in  the  two-element  case  that,  correct  to  first  order 
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effects  of  the  noise,  the  beamforming  technique  yields  an  estimate 
of  the  slope  that  has  the  same  variance  as  that  obtained  by 
making  independent  correlation  time  picks  in  each  waveform  and 
drawing  a  straight  line  between  the  resulting  times.  This 

observation  indicates  the  superiority  of  the  beamforming  tech- 

'•0  9 

nique  in  the  two-element  case  since  no  prior  assumption  about 
the  signal  is  necessary  to  carry  out  this  technique,  whereas  a 
perfect  replica  of  the  signal  is  used  in  the  correlation  time 
picks  involved  in  the  other  technique. 

The  two-element  case  has  been  presented  here  because 
the  equation  is  sufficiently  simple  that  it  is  easy  to  relate 
this  problem  to  that  of  the  correlation  time  pick.  The  N-element 
case  is  somewhat  more  complicated  (and  hence  less  transparent) 
and  is  considered  in  Appendix  B.  It  turns  out  that  just  as  in 
the  two-element  case,  assuming  large  signal-to-noise  ratios,  the 
beamforming  technique  and  the  plane-wave  fit  based  on  correlation 
time  picks  are  completely  equivalent.  Thus,  as  noted  above,  the 
beamforming  technique  appears  especially  valuable  in  the  case  of 
consistent,  but  unknown,  signal  waveshapes. 

The  derivation  in  Appendix  B  that  is  cited  in  the  preceding 
paragraph  assumes  a  perfectly  planar  wavefront  corrupted  by  additive 
noise.  The  complementary  derivation  also  appears  in  Appendix  B: 
no  noise,  but  random  small  perturbations  on  the  plane  wavefront. 

In  both  cases  the  beamforming  calculation  yields  the  same  result 
as  the  plane-wave  fit.  It  appears  that  the  same  result  would 
apply  (still  considering  only  first  order  effects)  in  the  combined 
.case  of  noise  and  travel-time  errors,  but  we  have  not  bothered 
to  demonstrate  this. 

3.5  EXAMPLES 

In  order  to  summarize  the  preceding  theoretical  results, 
we  shall  consider  some  specific  parameter  values.  If  we  assume 
a  uniform  N-element  array  of  total  aperture  200  km,  and  a  phase 
velocity  of  20  km/sec,  we  have  (see  Equation  3-3): 
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200 

m  =  4  sine  =  sine  =  7^-  sine  (3-22) 

2 

With  a  timing  error  variance  of  ot  we  would  have  (see  Equation 
3-11) 

^  in?  ~  i?  ? 

var  m  =  (rr-y)  var  sine  =  - 3 -  a  (3-23) 

iI_1  N(N2-1)  t 


If  we  assume  that  6  is  close  to  0, 
var  sine  =  var  6 
and,  therefore, 


var  e  = 


12 


N(N  -1) 


12  N-l 


100  NTN+1)  t 


(3-24) 


oQ  =  0. 35rtI-l/N(N+l)  ot  (3-25) 

To  be  more  specific,  the  radical  in  Equation  (3-25)  is  approxi¬ 
mately  0.28  (0.21)  for  N  =  11  (N  =  21).  Assuming  an  11-element 
array,  we  have 

og  «  0.1oT  (3-26) 

where  it  should  be  recalled  that  t  is  measured  in  seconds  and 
6  in  radians. 

If  the  t's  were  uniform  over  the  interval  0^1/40  sec, 
as  could  occur  as  a  result  of  unsynchronized  20/sec  sampling, 
we  would  have 

o  =  /3/3  •  1/40  =  0.0144 

T 

os  =  0.00144  radians  (3-27) 

y 

=  0.08° 

[With  the  21-element  array  this  number  would  be  0.06° . ] 

Using  Equation  (3-3)  we  can  calculate  what  noise  level 
would  contribute  the  same  amount  of  uncertainty  in  the  time  pick 
for  a  zero  crossing  time  pick.  Assuming  f()  =  lj  and  defining 
(S/N)  as  the  ratio  of  the  signal  energy  to  the  average  noise 
power,  we  have 
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(3-28) 


2 


(2*)  |(S/N)  2  tt  ( S/N ) 

The  ot  given  by  Equation  (3-27)  implies  a  (S/N)  given  by 

2 

(S/N)  =  3‘ ( ^  =  270  =  2k  dB  (3 

2  it 

As  another  example,  with  the  11-element  array.  Equa¬ 
tions  (3-28)  and  (3-26)  imply  that  a  10  dB  signal-to-noise 
ratio  leads  to 


o2  ~  0.005 

T 

os  -  o.n» 


(3-30) 


A  21-element  array  of  the  same  dimensions  would  imply 


(3-31) 
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SECTION  IV 


BEAMFORMING  AND  OTHER  METHODS  FOR  EPICENTER  LOCATION 

SIMULATION  RESULTS 

The  results  of  a  number  of  simulated  angular  measure¬ 
ments  using  the  center  elements  of  the  LASA  subarrays  are  described 
in  this  section.  As  described  in  Section  III,  the  analytical 
approach  toward  specifying  uncertainties  in  angular  measurements 
resulting  from  additive  noise  or  travel-time  anomalies  requires 
a  number  of  assumptions  that  may  not  be  realistic.  For  this  rea¬ 
son,  it  was  decided  that  a  better  understanding  of  the  effects 
of  these  kinds  of  interference  might  be  obtained  from  a  number  of 
computer-simulated  measurements.  To  perform  these  simulations, 
the  seismograms  from  a  magnitude-5.8  Kazakh  event  (as  recorded 
at  the  center  elements  of  the  twenty-one  LASA  subarrays)  were 
used  as  input  signals.  These  signals  were  chosen  for  two  reasons. 
First,  the  signal-to-noise  ratio  of  the  individual  seismograms 
was  high  (20  dB  or  greater);  by  adding  different  amounts  of  noise, 
it  was  therefore  possible  to  simulate  a  wide  range  of  signal- 
to-noise  ratios.  Moreover,  with  this  high  signal-to-noise  ratio, 
it  was  possible  to  remove  travel-time  anomalies  from  the  original 
records  prior  to  the  simulation  so  that  controlled  anomalies 
could  be  inserted  and  their  effects  studied.  Second,  this  col¬ 
lection  of  seismograms  is  typical  of  those  received  from  tele- 
seismic  ranges.  Since  the  records  differed  from  sensor  to  sensor, 
calculations  with  these  signals  are  expected  to  be  far  more 
realistic  than  those  which  could  be  made  assuming  identical  seis¬ 
mic  signals. 

To  simplify  calculations  these  seismograms  were  aligned, 
in  time,  so  that  they  appeared  to  have  come  from  an  epicenter 
due  north  of  the  array  at  a  range  corresponding  to  a  phase  velo¬ 
city  of  20  km/sec  (s80°).  No  attempt  was  made  to  compute  epicenter 
locations,  but  rather  the  effects  of  additive  noise  and  local 
travel-time  anomalies  on  arrival  angle  were  studied.  Further, 
computations  were  limited  to  the  case  of  a  constant  phase  velocity 
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(azimuth  scan),  since  the  results  obtained  in  this  manner  are 
directly  translatable  to  errors  in  range,  and  hence  to  errors 
in  epicenter  determination. 

This  discussion  begins  with  a  description  of  the  four 
methods  of  angular  measurement  that  have  been  studied:  plane- 
wave  fit,  analog  peak-signal  and  energy  beams,  and  a  DIMUS  beam. 
Following  this,  a  method  of  automatic  time-of-arrival  measure¬ 
ment  that  has  been  employed  as  an  input  for  all  of  the  angle 
measurement  system  mentioned  above  will  be  discussed.  The 
ground  rules  of  the  actual  computations  will  then  be  discussed, 
and  the  results  of  the  computations  presented.  Finally,  com¬ 
parisons,  where  possible,  will  be  made  with  the  theoretical  dis¬ 
cussion  of  Section  III. 

4 . 1  PLANE-WAVE  FIT 

The  following  procedure  was  used  to  determine  the  least- 
squares,  plane-wave  fit  to  the  arriving  signals.  First,  it  is 
assumed  that  the  general  angle  of  arrival  of  the  plane  wave  is 
known  to  within  a  small  error.  In  the  simulations,  the  correct 
arrival  azimuth,  due  north,  was  taken.  Next,  a  line  is  drawn 
through  the  center  element  of  the  array  (AO)  perpendicular  to 
the  assumed  direction  of  the  arrival  of  the  plane  wave.  Projec¬ 
tions  of  the  various  center  elements  of  the  LASA  clusters  on  this 
line  permit  the  determination  of  the  incremental  time  shifts  that 
must  be  applied  to  the  various  seismograms  in  "steering"  the 
array  in  the  vicinity  of  due  north.  Next,  arrival  times  of  a 
particular  feature  (in  this  case  a  zero  crossing  of  the  signal 
that  is  chosen  in  a  manner  that  will  be  described  below)  of  the 
signals  for  each  seismogram  are  compared  with  the  corresponding 
arrival  time  for  the  AO  seismogram.  Arrival  time  differences 
for  these  various  records  are  then  converted  to  distances,  using 
the  assumed  phase  velocity  of  20  km/sec  and  the  geometry  of  the 
projection  of  the  array.  A  least-squares  fit  of  these  distances 
to  a  line  passing  through  AO  gives  a  measure  of  the  difference 
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between  the  assumed  angle  of  arrival  and  that  which  is  estimated 
from  the  data.  Since  the  assumed  angle  of  arrival  is  the  correct 
one,  the  computation  gives  a  direct  measure  of  the  error  in 
angular  measurement  that  would  be  obtained  using  this  method. 

4 . 2  ANALOG  BEAM  -  PEAK  SIGNAL 

This  method  of  angular  measurement,  and  all  others  dis¬ 
cussed  in  this  report,  assumes  that  the  incoming  signal  can  be 
characterized  as  a  plane  wave.  With  this  assumption,  the  array 
is  steered  to  successive  angles,  and  the  maximum  value  of  the 
array  output  for  each  angle  is  taken  to  be  the  value  of  the  array 
pattern  for  that  angle  The  angle  of  arrival  of  the  signal  is 
defined  as  the  angle  for  which  this  pattern  Assumes  its  maximum  value. 

A  more  detailed  description  of  this  calculation  follows. 
First,  the  seismograms  from  each  element  of  the  array  are  nor¬ 
malized  so  that  the  peak  of  the  absolute  value  of  the  signals 
is  identical  for  all  sensors.  The  various  signals  are  then  ad¬ 
vanced  or  delayed  in  time.  These  delays  are  computed  from  the 
projection  of  the  location  of  the  sensors  on  a  line  through  the 
center  of  the  array,  the  steering  angle,  and  the  assumed  phase 
velocity  of  20  km/sec.  The  seismograms  are  then  summed,  and 
the  sum  is  examined  in  a  time  window  of  two  seconds  centered  on 
the  arrival  time  of  the  signal  at  the  center  element  of  the 
array.  The  peak  of  the  absolute  value  of  the  sum  seismogram 
in  this  time  window  is  taken  to  be  the  array  output.  The  choice 
of  time  window  was  made  to  include  the  coherent  portions  of  the 
seismograms,  and  to  allow  for  some  advance  and  retardation  of 
signals  from  sensors  other  than  the  one  at  the  center  of  the 
array . 

For  all  pattern  calculations  described  in  this  section, 
the  patterns  were  developed  in  0.1°  steps  for  the  region  of  *3° 
about  the  assumed  angle  of  arrival.  In  an  attempt  to  minimize 
fine  structure  in  the  main  lobe  of  the  array  patterns  obtained, 
local  averaging  of  the  pattern  was  performed.  Each  point  of  the 
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patterns  presented  represents  the  average  of  the  five  adjacent 
points  of  which  it  is  the  center.  No  analytical  optimization  of 
this  smoothing  was  attempted,  and  in  fact,  the  experimental  results 
suggest  that  further  smoothing  might  produce  better  results. 

4.3  ANALOG  BEAM  -  ENERGY 

This  computation  of  the  angle  of  arrival  of  signals 
corresponds  exactly  to  that  just  described,  except  that  the  beam 
output  is  taken  to  be  the  energy  contained  in  the  time  window 
rather  than  the  peak  signal  in  this  interval.  This  is  essentially 
the  beamforming  calculation  that  is  considered  in  Section  III. 

This  method  was  expected  to  perform  somewhat  better  than  the  peak 
or  plane-wave  fit  in  the  presence  of  low  signal-to-noise  ratios. 

4.4  DIMUS  BEAM 

With  a  DIMUS  processor,  the  outputs  of  the  sensors  of 
the  array  are  hardlimited  before  any  array  processing  takes  place. 
Thus  the  array  output  consists  of  the  summation  (properly  delayed 
in  time)  of  the  hardlimited  outputs  of  the  various  sensors  of  the 
array.  The  peak  value  definition  of  a  pattern  would  not  be  satis¬ 
factory  in  this  case,  since  the  peak  value  of  this  summed  seis¬ 
mogram  in  the  time  window  described  above,  is  constant  over  the 
range  of  angles  which  fall  between  the  3  dB  points  of  the  analog 
(peak)  beam  pattern.  This  is  because  the  relative  time  delays 
that  are  required  across  the  array  to  cause  the  conventional  array 
pattern  to  fall  to  the  3  dB  level  correspond  roughly  to  a  relative 
time  offset  of  one-half  period  between  sensors  located  at  the 
extremes  of  the  aperture  of  the  array.  With  hardlimited  signals, 
this  is  precisely  the  minimum  time  shift  which  can  cause  any 
change  in  the  array  output,  since  the  hardlimited  signals  remain 
constant  over  each  half  cycle  of  the  signals. 

Because  of  this  property  we  have  devised  an  alternate 
calculation  for  the  DIMUS  array  measurements  of  angle,  in  which 
the  conventional  hardlimited  outputs  of  the  sensors  are  replaced 
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with  synthetic  seismograms  which  are  zero  except  in  the  vicinity 
of  the  zero  crossings  of  the  original  seismograms.  At  those 
sample  points  at  which  the  original  seismogram  possesses  a  zero 
crossing,  that  sample  point,  and  the  one  immediately  succeeding, 
are  set  to  unity  with  a  sign  which  corresponds  to  the  direction 
of  the  zero  crossing  of  the  original  seismogram.  These  synthetic 
seismograms  are  then  delayed,  as  in  the  case  of  steering  an 
analog  beam,  and  summed.  This  summed  seismogram  is  then  squared, 
and  the  integral,  with  respect  to  time,  is  taken  over  the  two- 
second  interval  centered  around  the  onset  of  the  arrival  time  of 
the  seismic  signal. 

Since  the  total  information  of  the  hardlimited  inputs 
to  the  DIMUS  must  be  contained  in  the  zero  crossings  of  the 
signals,  this  form  of  processing  gives  large  outputs  when  zero 
crossings  of  the  same  sign  are  coincident  in  the  synthetic 
seismograms  (n  vs.  n,  where  n  is  the  number  of  (coincident) 
zero  crossings),  and  also  gives  a  premium  for  alignment  of 
several  zero  crossings  in  sequence  (a  factor  of  m,  where  m  is 
the  number  of  successive  zero  crossings  within  the  time  window) . 

k . 5  AUTOMATIC  DETECTION  OP  ARRIVAL  TIME 

The  various  methods  of  angular  determination  described 
in  the  preceding  sections  all  require  a  knowledge  of  the  onset 
time  of  the  signal  at  the  array.  Although  this  arrival  time  is 
known  a  priori  from  the  construction  of  the  simulation  tests, 
it  was  believed  that  the  simulations  would  be  more  realistic  if 
they  included,  as  an  integral  part  of  the  calculations,  the 
automatic  detection  of  arrival  time  as  measured  from  the  seismo¬ 
grams  themselves.  The  implementation  of  this  automatic  detector 
was  carried  out  with  the  output  of  the  unphased  sum  of  a  DIMUS 
beam  as  discussed  in  a  previous  report  [  7  ].  It  has  earlier  been 
shown  that  with  DIMUS  processing,  it  is  very  likely  that  the  first 
motion  of  seismic  signals  will  be  emphasized,  and  in  fact  nearly 
always  will  result  in  DIMUS  outputs  that  saturate  during  the  first 
few  half  cycles  of  the  seismic  record  [8].  Thus,  in  the  case  of 
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the  simulation  tests,  it  is  expected  that  the  portion  of  the  DIMUS 
output  corresponding  to  first  motion  will  be  nearly  equal  to  21 
(the  number  of  seismometers  in  the  array).  Hence,  if  the  magnitude 
of  the  DIMUS  output  exceeds  a  threshold  set  at  18,  the  following 
zero  crossing  should  correspond  to  that  which  we  have  used  as 
a  definition  of  arrival  time.  This  zero  crossing  on  the  sum 
seismogram  is  then  used  as  a  reference  point,  and  the  individual 
seismograms  are  searched  in  the  vicinity  of  this  time  for  zero 
crossings  of  the  same  polarity  as  that  for  the  array  output. 

These  zero  crossings  are  then  identified  as  the  arrival  times  of 
the  signals  at  the  various  sensors . 

Of  course,  the  output  of  the  DIMUS  will  be  a  function 
of  the  direction  in  which  the  array  is  steered;  however,  as 
previously  stated,  the  peak  signal  output  from  the  array  is 
expected  to  be  only  a  weak  function  of  angle  in  the  vicinity  of 
the  actual  arrival  angle  of  the  wavefront.  To  test  the  validity 
of  this  assumption,  a  linear  array  (11  elements)  with 
a  200  km  aperture  was  synthesized  using  eleven  of  the  seismograms 
used  for  the  rest  of  the  simulation  tests,  and  the  array  was 
scanned  in  three-degree  steps  with  a  signal-to-noise  ratio  of 
3  dB.  Results  of  this  simulation  are  shown  in  Figure  4.1.  From 
the  figure,  it  can  be  seen  that,  provided  the  DIMUS  beam  is 
steered  within  a  small  region  (several  degrees)  of  the  true 
arrival  angle,  the  detector  should  respond  properly.  Hence, 
the  unphased  DIMUS  detector  appears  to  work  well,  and,  since 
the  automatic  detection  of  arrival  times  makes  for  a  more  realis¬ 
tic  simulation,  it  was  employed  in  all  of  the  calculations 
presented  here. 

Before  proceeding,  it  should  be  noted  that  the  DIMUS 
detector  did  not  always  identify  the  proper  zero  crossing  in 
the  simulations.  On  occasion  it  picked  either  the  actual  onset 
of  first  motion  of  the  signals  or  the  second  zero  crossing 
after  this  time.  Examination  of  results  obtained  using  only 
the  proper  zero  crossing,  as  opposed  to  using  all  results. 
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FIGURE  4.1 

DIMUS  ARRAY  OUTPUTS  FOR  SELECTED  AZIMUTHS 
SIGNAL  TO  NOISE  «3dB 
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showed  no  significant  difference  in  the  angular  measurements 
to  be  discussed  in  later  sections.  As  a  result,  no  further 
attention  was  paid  to  this  phenomenon. 

4.6  METHOD  OF  ADDING  NOISE  AND  ANOMALIES 

The  noise  samples  added  to  the  various  seismograms 
to  test  the  effects  of  signal-to-noise  ratio  on  angular  accuracy 
were  developed  in  the  following  manner.  A  random  number  generator 
was  used  to  develop  twenty-one  uncorrelated  series  of  numbers  for 
each  simulation  trial.  These  series  of  numbers  were  taken  to  be 
samples  of  noise  taken  0.05  seconds  apart.  To  tailor  the  spectrum 
to  these  noise  signals  to  represent  an  approximation  of  seismic 
noise,  each  series  was  passed  through  a  recursive  filter,  derived 
from  a  three-pole  Butterworth  filter,  with  3  dB  points  at  0.5 
and  1.5  Hz,  The  twenty-one  independent  noise  records  thus 

obtained  were  scaled  by  an  appropriate  constant  and  added  to  the 
twenty-one  input  seismograms.  As  stated  earlier,  these  seismograms 
are  normalized  so  that  the  peak  value  of  each,  in  the  time  window 
from  first  motion  to  one  second  later,  is  equal  to  unity.  The 
noise  scale  factor  was  obtained  for  each  signal-to-noise  ratio 
by  assuming  that  each  of  the  original  seismograms  had  the  same 
signal  energy  as  one  period  of  a  1  Hz  sinewave  of  unity  peak  value 
for  two  trials  —  one  at  10  dB,  and  one  at  3  dB  signal-to-noise 
ratio.  Successive  trials  were  then  made  with  different  sets  of 
noise  samples.  Twelve  such  pairs  of  trials  were  included  in  the 
simulations  reported  here. 

For  the  simulations  used  to  study  the  effects  of  local 
travel-time  anomalies,  the  following  procedure  was  used.  For  each 
trial,  a  set  of  twenty  independent  random  numbers  was  generated 
from  a  uniform  distribution  defined  between  in/20  sec,  where  n 
is  the  maximum  error  in  sample  points.  The  resulting  numbers 
were  rounded  off  to  the  nearest  number  of  sample  points  (between 
in)  and  applied  as  delays  to  each  of  the  seismograms,  except  AO. 
Nine  simulations  each  were  made  for  n=l  and  n=2.  It  should  be 
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noted  that  the  resulting  (discrete)  anomaly  distributions  are  not 
quite  uniform.  (The  probabilities  associated  with  in  are  too  small 
by  a  factor  of  1/2.)  The  standard  deviations  of  the  anomalies 
generated  as  described  were  0.035  and  0.07  seconds,  respectively. 

4.7  RESULTS  WITH  NOISE  ALONE 

Results  of  the  simulations  with  noise  alone  can  be  illus¬ 
trated  by  Figures  4.2,  4.3  and  4.4.  Figure  4.2  gives  representative 
seismograms  at  element  Bl.  Figure  4.3  shows  the  outputs  of  the 
analog  and  DIMUS  beamformers  for  one  trial  with  a  signal-to-noise 
ratio  equal  to  10  dB ,  Also  illustrated  in  this  figure  are  the 
analog  array  peak  and  energy  patterns,  as  well  as  the  DIMUS  array 
pattern.  On  these  patterns,  a  +  indicates  the  peak  of  the  pattern 
and  a  A,  the  angular  location  of  the  plane-wave  fit.  Figure  4.4 
is  the  same  presentation  for  a  signal-to-noise  ratio  of  3  dB . 

Figure  4.5  presents  histograms  of  the  number  of  trials  with  a 
given  error  in  angular  measurement  for  each  of  the  four  computation 
schemes  employed:  plane-wave  fit,  analog  peak  and  energy  patterns, 
and  the  DIMUS  pattern.  Both  the  3  and  10  dB  cases  are  included. 

From  the  figure  it  can  be  stated  that  the  distribution  of  errors 
appears  nearly  uniform  for  the  lower  signal-to-noise  ratio,  and 
approaches  a  more  peaked  shape  for  the  higher  signal-to-noise  ratio. 

Table  I  summarizes  these  results  in  terms  of  the  means 
and  standard  deviations  of  the  angle  estimates  without  regard  to 
a  priori  knowledge  of  the  actual  arrival  angle.  Examination  of 
the  table  permits  several  observations.  First,  there  appears 
to  be  a  bias  in  the  estimate  of  the  arrival  angle,  which  shows 
a  general  decrease  with  increasing  signal-to-noise  ratio.  Second, 
the  standard  deviation  of  the  error  in  measurement  of  arrival  angle 
(from  the  mean  of  the  data)  generally  decreases  with  increasing 
signal-to-noise  ratio.  Finally,  the  errors  in  measurement  obtained 
with  a  signal-to-noise  ratio  of  10  dB  are  within  acceptable  limits 
for  azimuth  determination  in  an  operational  system. 
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FIGURE  4.2 

SIGNAL  SEISMOMETER  OUTPUTS  FOR  NINE  SIMULATION  TRIALS 
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FIGURE  4.3 

REPRESENTATIVE  RESULTS  FOR  ONE  TRIAL  WITH  NOISE  ADDED 

SIGNAL  TO  NOISE  *  10  dB 
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FIGURE  4.4 

REPRESENTATIVE  RESULTS  FOR  ONE  TRIAL  WITH  NOISE  ADDED 

SIGNAL  TO  NOISE  =  3  dB 
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TABLE  I  -  SAMPLE  MEANS  AND  STANDARD  DEVIATIONS  VS.  NOISE 
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TABLE  II  -  SAMPLE  MEANS  AND  STANDARD  DEVIATIONS  VS.  TIME  ANOMALIES 
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4.8  RESULTS  WITH  TIME  ANOMALIES  ALONE 

Results  obtained  from  simulations  with  travel-time 
anomalies  are  illustrated  by  Figures  4.6  and  4.7.  For  the  case 
of  a  one-sample  point  maximum  anomaly  Figure  4.6  gives  the 
summed  seismogram  and  the  DIMUS  sum  (both  were  obtained  with 
the  array  steered  at  the  known  arrival  angle),  the  analog 
pattern-peak,  the  analog  pattern-energy,  and  the  DIMUS  pattern. 
Shown  on  the  plots  of  the  array  patterns  are  the  maximum  value 
of  the  pattern,  indicated  by  a  +,  and  the  angle  of  the  plane 
wave  fit,  indicated  by  a  A.  Figure  4.7  gives  the  same  infor¬ 
mation  for  the  case  of  a  two-sample  point  maximum  anomaly. 

Figure  4.8  presents  histograms  of  the  number  of  trials 
with  given  angular  errors  for  each  of  the  four  methods  and  for 
both  the  cases  of  one-  and  two-sample  point  maximum  anomaly. 

These  results  are  summarized  in  Table  II.  Examination  of 
Table  II  shows  that  in  general  the  mean  error  and  the  standard 
deviation  of  this  error  increase  with  the  size  of  the  anomalies. 
As  in  the  noise-alone  case,  angular  errors  that  would  be  adequate 
for  systems  applications  appear  feasible  in  the  less  severe  of 
the  two  cases  considered. 

4.9  CONCLUSIONS 

Examination  of  the  histograms  of  errors  for  all  simu¬ 
lation  trials  made  (Figures  4.5  and  4.8)  indicates  insufficient 
information  to  make  more  than  the  grossest  of  estimates  of  the 
distribution  of  these  'errors.  Taken  together  with  the  fact  that 
the  theoretical  treatment  of  Section  III  dealt  with  continuous, 
rather  than  discrete  times,  it  is  possible  only  to  make  quali¬ 
tative  comparisons  of  predicted  and  simulated  results.  With 
this  in  mind,  compare  the  standard  deviations  in  angle  obtained 
for  a  10  dB  signal-to-noise  ratio  with  the  value  of  0.3°  pre¬ 
dicted  by  theory.  The  spread  in  simulation  results  lies  between 
0.18°  and  Q.33°3with  the  smallest  standard  deviation  belonging 
to  a  population  with  a  substantial  mean  error.  The  agreement  is 
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PEAK  ANOMALY  =  I  SAMPLE  POINT 
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NUMBER  OF  TRIALS  WITH  INDICATED  ANGLE  ERRORS 
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NUMBER  OF  TRIALS  WITH  INDICATED  ANGLE  ERRORS 


FIGURE  4.8  CONT. 
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as  good  as  could  be  expected.  Examination  of  Table  I  shows 
that  the  standard  deviations  of  error  increase  with  decreasing 
signal-to-noise  ratio.  The  mean  of  the  ratios  of  the  3  dB  to 
10  dB  standard  deviations  is  2.5,  which  compares  quite  well  with 
the  theoretically  predicted  ratio  of  2.2.  If  the  DIMUS  beam 
result  is  removed  from  the  set  (there  is  no  reason  to  expect 
the  peculiar  calculations  of  this  beam  to  agree  with  the  theore¬ 
tical  treatment  of  Section  III),  the  mean  ratio  is  exactly 
2.2  --  better  than  should  be  expected. 

Results  of  the  anomaly  simulations  give  standard  devi¬ 
ations  of  errors  which  are,  in  general,  higher  than  theoretical 
predictions  by  a  factor  of  nearly  two  for  both  the  one  and  two 
sample  point  peak  errors  (see  Table  II).  What  is  more  important, 
there  are  biases  in  the  data  samples  producing  substantial 
(0.2°)  mean  errors  in  angular  determination.  It  is  hypothesized 
that  both  of  the  above  effects  are  due  to  the  poor  estimates 
of  the  error  distributions  available  from  the  simulation  trials. 
What  is  perhaps  most  important,  is  that  it  is  unlikely  that 
travel-time  anomalies  can  be  reduced  below  the  one  sample  point 
level.  Thus,  while  the  results  obtained  suggest  that  for  large 
enough  events,  noise  will  not  substantially  limit  the  accuracy 
of  event  location,  we  may  always  be  left  with  a  standard  devia¬ 
tion  of  angular  error  on  the  order  of  0.25°.  Since  this  result 
is  in  such  large  disagreement  with  the  theoretically  predicted 
one,  further  study  is  indicated.  Before  resolution  of  this 
problem,  however,  the  simulation  results  may  be  translated  into 
geographical  terms.  This  would  imply  a  standard  deviation  of 
some  25  kilometers  in  one  component  of  the  determination  of  an 
epicenter  from  a  single  large  array. 
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SECTION  V 


CODA-CORRELATION  STATISTIC  WITH  DIMUS  (HARDLIMITED)  WAVEFORMS 

In  a  previous  report  [9  ]>  data  were  presented  that 
indicated  the  average  paired  correlation  coefficient  could  be 
a  useful  discriminant  for  identifying  earthquakes  with  depths 
between  10  and  40  km.  These  calculations,  which  were  based  on 
the  outputs  of  seismometers  distributed  over  thousands  of 
kilometers,  indicated  significantly  more  coda  correlation  for 
earthquakes  in  the  10-40  km  region  than  for  other  earthquakes 
or  for  surface-focus  events. 

We  have  been  studying  the  possibility  of  using  DIMUS 
(hardlimited)  seismograms  for  two  other  discriminants,  com¬ 
plexity  and  spectral  ratio  [10],  and  for  location  [ll].  Along 
with  these  studies,  it  is  of  interest  to  examine  how  much  the 
coda-correlation  discriminant  would  be  degraded  by  the  use  of 
DIMUS  waveforms.  Results  of  repeating  these  calculations  with 
DIMUS  waveforms,  obtained  from  all  but  one  of  the  same  set  of 
42  events,  are  summarized  in  this  section.  It  turns  out  that 
the  discriminant  still  appears  useful  for  separating  "moderate 
depth"  earthquakes  from  other  events,  but  the  separation  is 
not  as  clear  as  it  was  in  the  case  of  analog  waveforms. 

The  calculation  of  the  average  paired  correlation 
coefficient  has  been  described  in  detail  elsewhere  [12].  For 
the  DIMUS  calculations,  since  all  waveforms  have  a  unity  rms 
value  in  all  time  intervals,  an  alternate  (but  equivalent)  imple¬ 
mentation  of  the  calculation  turns  out  to  be  much  simpler.  This 
implementation  is  based  on  a  simple  relation  between  the  average 
paired  correlation  coefficient  and  the  energy  in  the  sum  of  the 
normalized  waveforms. 

Figure  5.1  is  a  scatter  plot  of  the  coda  correlation 
coefficients  calculated  from  the  analog  waveform,  p^,  and  those 
calculated  from  the  DIMUS  waveform,  p^.  Different  symbols  are 
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used  for  the  three  different  classes  of  events:  earthquakes 

with  depths  between  10  and  40  km  (23  events),  other  earthquakes 

(6  events),  and  surface-focus  events  (12  events).  Two  lines 

are  drawn  to  illustrate  the  separation  of  the  points  obtained 

for  the  moderate  depth  earthquakes  from  the  others.  With  the 

exception  of  one  surface-focus  event  (with  =  0.2),  it  is 

clearly  easy  to  separate  the  moderate  depth  earthquakes  from 

the  other  events  with  the  horizontal  line  at  p  =  0.14.  In  the 

c 

case  of  the  DIMUS  discriminant,  a  vertical  line  at  p,  =  0.05 
separates  most,  but  not  all  of  the  events,  and  there  are  several 
points  very  close  to  the  line. 

It  is  clear  that  the  DIMUS  coda-correlation  statistic 
is  inferior  to  the  statistic  calculated  from  the  analog  wave¬ 
forms.  However,  it  still  provides  a  fairly  good  separation 
of  moderate  depth  earthquakes  from  other  events.  We  expect, 
in  the  near  future,  to  be  processing  additional  events  with  the 
analog  coda-correlation  statistic.  For  this  reason  we  shall 
postpone  any  further  comments  or  conclusions  on  the  performance 
of  the  DIMUS  coda-correlation  statistic  until  these  additional 
calculations  have  been  performed. 
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SECTION  VI 


EXTENDED  AUTOMATIC  pP  TEST 

The  previously  reported  work  [13]  on  an  extended  automatic 
pP  test  has  been  continued  using  arrays  of  continental  dimensions. 

It  was  anticipated  that  use  of  continental-size  arrays  would  yield 
good  signal-to-coda  enhancement  in  the  phased  sum  seismogram  due  to 
coda  decorrelation  over  the  large  spatial  separations  of  the  seis¬ 
mometers,  and  due  to  the  increased  moveout  of  the  pP  phase  across 
the  array.  The  key  to  the  success  or  failure  of  the  test  was  the 
pP  alignment  procedure,  which  is  based  on  an  assumed  depth  of  the 
Moho  and  assumed  average  velocities  above  and  below  the  Moho .  The 
P-pP  delay  times  for  alignment  were  calculated  using  a  simple  geo¬ 
metric  model  of  surface  reflection  described  below.  It  was  pre¬ 
viously  indicated  and  should  be  emphasized  that  errors  of  10  to 
20  km  in  depth  determination  of  deep-focus  events  are  not  important 
for  discrimination  purposes. 

To  date,  the  extended  test  has  been  applied  to  eight  events 
using  arrays  of  continental  dimensions  with  quite  good  results. 
Information  concerning  these  events  is  summarized  in  Table  III. 

In  several  of  these  events,  the  test  detected  sP  as  well  as  pP, 
which  proved  to  be  a  valuable  aid  in  depth  determination. 

6.1  DESCRIPTION  OP  TEST 

As  described  earlier,  the  automatic  pP  test  proceeds  by 
calculating  for  each  of  several  "test  depths"  a  delayed  sum,  in 
which  each  seismogram  is  aligned  in  such  a  way  that  the  pP-phases 
would  add  coherently  if  the  test  depth  were  correct.  Then  the  en¬ 
ergy  in  the  one-second  time  window  that  would  enclose  the  summed 
pP-phase  is  calculated  and  divided  by  a  measure  of  the  surrounding 
energy.  Finally,  this  energy  ratio  statistic  is  multiplied  by  an 
average  paired  correlation  coefficient  and  the  test  depth  yielding 
the  largest  value  is  selected. 
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TABLE  III  -  EVENTS  USED  FOR  AUTOMATIC  pP  TEST 
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Three  changes  have  been  made  in  converting  the  test  from 
the  10-40  km  region  to  deeper  values.  The  first  change,  which  is 
probably  not  very  important,  is  that  the  denominator  of  the  maximal- 
ratio  weights  used  in  combining  seismograms  is  now  calculated  from 
a  50-second  time  window  following  the  P  wave,  rather  than  the  10- 
second  window  used  previously.  The  second  change  is  in  the  calcu¬ 
lation  of  the  energy  ratio.  The  energy  in  the  one-second  test- 
depth  window  is  now  compared  to  the  average  energy  in  a  five-second 
window  preceding  the  test-depth  window.  The  third  change  relates 
to  the  alignment  procedure,  and  is  crucial  to  the  success  of  the 
extended  test. 

The  P-pP  delay  times  are  calculated  from  the  JB  travel-time 
tables  using  the  geometry  illustrated  in  Figure  6.1.  Point  E  is 
the  hypocenter  of  the  event;  h  is  the  depth  of  the  event.  The  pP 
phase  is  reflected  at  point  0  on  the  surface  of  the  earth,  and 
passes  through  point  R  at  depth  h  along  its  path  to  the  seismometer. 
The  incident  angle  1q  is  calculated  from  travel-time  tables  and  is 
almost  invariant  with  depths  down  to  200  km.  The  Moho  is  indicated 
in  Figure  6.1  at  a  depth  of  h  ,  and  all  of  the  following  calculations 
are  made  for  h  greater  than  h  . 

The  alignment  times  are  computed  as  follows.  The  P-phase 
travel  time  to  the  seismometer  is  calculated  from  the  JB  tables, 
interpolating  linearly  with  range  and  depth.  The  pP-phase  travel 
time  to  the  seismometer  is  found  in  two  steps  by  adding  the  times 
of  travel  from  E  to  0  to  R,  say  t^,  and  from  R  to  the  receiver. 

Let  v.  and  v.  denote  the  average  velocities  above  and  below  the 
Moho  discontinuity  respectively,  then 


h 

t,  =  2  [—  + 
1  va 


(h-hj 


v. 


-]/cos(1q)  seconds 


(6.1) 


Possible  refraction  effects  at  velocity  discontinuities  are  neglected. 
The  travel  time  from  R  to  the  receiver  is  from  JB  tables  as  above, 
only  the  range  is  now  the  range  of  E  minus  ER,  where 

ER  =  2h  tan(iQ)  km  =  2h  tan( iQ ) /111 . 195  degrees  (6.2) 
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FIGURE  6.1 

GEOMETRY  OF  SURFACE  REFLECTION 
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The  parameter  values  =  6.34  km/sec,  =  8.00  km/sec,  and 

h  =  33.0  km  have  been  used  in  all  calculations  to  date.  Moderate 
m 

errors  in  the  assumed  parameter  values  will  not  seriously  affect 
the  results  of  the  test.  For  example,  if  the  Moho  is  misplaced 
by  30  km,  or  if  the  assumed  average  velocities  are  in  error  by 
10%,  then  the  two  extreme  phasing  times  will  be  in  error  by  less 
than  0.10  sec,  for  an  event  with  100  km  depth,  with  epicenter 
approximately  40°  from  an  array  with  a  2500  km  aperture.  This  is 
due  to  the  fact  that  the  P-wave  incident  angles  at  the  receivers 
are  essentially  invariant  with  depth.  Although  these  parameter 
variations  would  not  lead  to  serious  phasing  errors,  they  could 
cause  the  test  statistic  to  peak  at  the  wrong  depth,  which,  as 
noted  above,  would  not  be  critical  for  discrimination  purposes. 

An  inherent  property  of  the  automatic  pP  test  is  that 
it  may  also  detect  the  sP  phase  arrival.  The  sP  phase  delay 
relative  to  the  P-phase  is  a  constant  factor  times  the  pP  phase 
delay.  Since  the  pP  phase  delays  vary  linearly  with  depth,  when 
the  automatic  test  calculates  the  test  statistic  at  a  depth 
approximately  one-and-one-half  times  the  depth  of  the  event,  the 
sP  phase  signals  will  add  coherently,  and  could,  if  the  signals 
were  strong  enough,  cause  the  test  statistic  to  peak. 

For  some  events,  poor  signal-to-noise  ratios  necessitated 
filtering  of  the  data.  The  test  results  for  two  events  were 
vastly  improved  after  filtering,  which  was  done  with  a  recursive, 
bandpass  digital  filter  derived  from  a  low-pass  three-pole  Butter- 
worth  analog  filter.  The  3  dB  points  of  the  passband  were  0.6  Hz 
and  3.25  Hz. 

6 . 2  EXPERIMENTAL  RESULTS 

A  measure  of  the  effectivenss  of  the  automatic  depth 
test  might  be  a  comparison  with  the  C&GS  reported  depths.  C&GS 
depths,  though .  seem  not  to  be  a  reliable  standard,  and  whenever 
possible,  we  have  referred  to  the  original  seismograms  to  resolve 
discrepancies.  Only  two  of  the  eight  events  tested  agreed  with 
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the  C&GS  depth  determinations.  Two  other  events  suggest  that 
C&GS  mistakenly  identified  sP  as  pP.  The  depth  results  of  three 
other  events  conflict  with  the  C&GS  reported  depths  rather 
severely.  A  presentation  of  seismograms  with  phase  identifica¬ 
tions  is  made  in  support  of  the  automatic  test  for  one  of  these 
three  events,  and  for  the  two  events  where  we  believe  sP  was 
mistaken  for  pP .  The  final  event  was  listed  by  C&GS  as  33  km, 
indicating  that  they  were  unable  to  make  a  depth  estimate;  the 
automatic  test,  however,  did  show  a  peak  in  the  test  statistic, 
which  is  again  supported  by  illustrative  seismograms. 

The  two  events  which  demonstrated  good  agreement  between 
the  automatic  test  depths  and  the  C&GS  depths  both  had  large 
signal-to-noise  ratios.  The  21  March  1962  event  had  an  average 
single  station  signal-to-noise  ratio  of  30,  and  the  15  July  1963 
event  had  an  average  single  station  signal-to-noise  ratio  of  212. 

In  both  cases,  the  phases  were  easily  identifiable  by  eye  on  all 
seismograms.  The  test  results  are  shown  in  Figure  6 . 2D  and 
Figure  6.2B,  respectively. 

The  two  events  for  which  we  believe  sP  was  mistaken  as 
pP  by  C&GS  were  12  September  1962  and  19  July  1962.  The  test 
statistic  of  Figure  6.2A  shows  the  pP  phase  clearly  at  101  km 
and  the  sP  phase  following  at  153  km  for  the  12  September  1962 
event.  This  delay  corresponds  well  to  relative  transverse  and 
longitudinal  wave  velocities.  An  illustrative  set  of  seismograms 
is  presented  in  Figure  6.4,  with  the  pP ,  sP  and  PcP  phases 
identified  on  each  at  locations  corresponding  to  a  depth  of 
approximately  101  km.  The  19  July  1962  event  produced  only 
one  rather  strong  peak  when  the  automatic  test  was  performed  with 
unfiltered  seismograms  (see  Figure  6.2C).  This  peak  occurred 
at  a  depth  of  147  km,  13  km  from  the  C&GS  listed  depth,  and 
would  have  been  accepted  as  the  correct  depth  except  that  the 
test  was  repeated  with  filtered  seismograms.  After  filtering, 
another  clear  peak  emerged  at  97  km,  as  can  be  seen  in  Figure  6 . 3A , 
indicating  that  the  peak  at  147  km  was  in  fact  the  sP  phase. 
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A  review  of  the  seismograms  for  this  event  indicates  that  97  km 
is  the  correct  depth  for  this  event.  Seismograms  for  this  event 
are  presented  in  Figure  6.5  with  the  pP  and  sP  arrival  times  indi 
cated  for  a  97  km  depth  event 

The  three  events  for  which  the  automatic  test  disagrees 
with  C&GS  estimates  occurred  on  25  August  1963,  10  January  1963, 
and  6  October  1962.  The  first  of  these,  the  25  August  1963  event 
was  studied  in  an  earlier  report  [l4]  and  seemed  to  be  an  event 
with  depth  between  10  and  40  km,  although  C&GS  has  it  listed  at 
13^  km.  When  the  extended  test  was  performed  with  this  event, 
there  were  no  peaks  in  the  test  statistic  in  the  33  to  200  km 
depth  range.  The  maximum  value  of  the  test  statistic  was  0.35* 

An  isolated  peak  of  magnitude  less  than  1.0  in  the  test  statistic 
is  not  considered  meaningful  by  itself;  hence  we  believe  that 
this  event  was  not  in  the  33  to  200  km  depth  region.  No  plot 
is  presented  for  this  event. 

The  10  January  1963  event  was  more  interesting.  The 
automatic  test  showed  a  depth  of  55  km  (see  Figure  6.2E),  while 
the  C&GS  depth  is  125  km.  Three  illustrative  seismograms  are 
presented  in  Figure  6.6,  showing  pP,  sP  and,  on  one  seismogram, 
PcP  phases  for  an  event  of  55  km  depth.  The  seismograms  are 
rather  noisy  and  yet  the  pP  phases  corresponding  to  a  55  km 
depth  event  stand  out  quite  well. 

The  6  October  1962  event  was  a  low  magnitude  event  (the 
actual  magnitude  is  not  listed  by  C&GS) .  The  pP  phase  could 
not  be  identified  by  C&GS  at  any  station.  The  C&GS  depth  of 
149  km  comes  from  an  unrestrained  epicenter  calculation,  which, 
however,  was  supported  by  S-phase  arrivals  at  several  stations 
at  close  proximity  to  the  event  location  [15].  The  seismograms 
for  this  event  are  very  poor,  and  the  automatic  test  yielded 
no  peaks  in  the  test  statistic  when  performed  with  unfiltered 
data.  When  performed  with  filtered  data,  the  automatic  test 
showed  a  vast  improvement.  As  shown  in  Figure  6.3B,  a  clear 
peak  occurred  at  77  km,  and  a  smaller  peak  occurred  at 
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109  km,  suggesting  pP  and  sP,  respectively.  We  have  not  yet  had 
the  opportunity  to  confirm  this  result  on  the  filtered  seismograms. 

The  final  event,  6  January  1964,  is  a  magnitude  5.6 
earthquake,  which  overloaded  several  of  the  stations  used  in  the 
automatic  depth  test.  In  addition,  the  seismograms  (see  Figure  6.7) 
show  very  strong  codas,  which  make  it  easy  to  misidentify  a  phase 
arrival.  As  a  result,  C&GS  made  no  depth  estimate  for  this  event. 
The  automatic  test  yielded  a  clear  peak  at  155  km,  which  is  diffi¬ 
cult  to  verify  or  refute  from  the  seismograms.  Overloading  is 
present  on  three  of  the  nine  stations,  and  this  badly  degrades  the 
quality  of  those  signals.  If  there  were  a  phase  arrival  at  a  point 
of  overloading,  the  signal  distortion  would  cause  a  low  average 
paired  correlation,  and  would  probably  cause  the  automatic  test  to 
miss  that  phase.  In  support  of  the  automatic  depth  test  peak  at 
155  km,  it  is  observed  that  the  average  paired  correlation  was 
0.35  and  the  energy  in  the  test  depth  window  compared  to  the 
average  energy  in  the  coda  window  was  4.7.  An  extraneous  peak 
of  this  nature  is  highly  unlikely  in  the  signal  region  following 
the  P-phase  by  more  than  ten  seconds.  Since  the  only  phase  arrivals 
which  could  possibly  add  coherently  in  the  automatic  test  are  pP 
and  sP,  it  seems  likely  that  this  event  occurred  at  a  depth  of  at 
least  100  km.  The  results  of  the  test  are  shown  in  Figure  6.2F. 

6.3  SUMMARY  AND  CONCLUSIONS 

The  extended  automatic  pP  test  has  been  applied  to  eight 
events  using  continental-size  arrays.  The  experimental  results 
were  good  for  all  eight  events.  In  three  events  the  test  detected 
the  sP  phase  arrivals  as  well  as  the  pP  phase  arrivals.  On  two 
of  these  three  events  it  appears  that  C&GS  also  detected  the  sP 
phase  arrivals  and  mistook  them  for  pP.  Since  sP  delay  times  are 
the  same  as  the  pP  delay  times  within  a  constant  factor,  both 
phases  can  be  added  coherently  by  the  automatic  test.  The  appear¬ 
ance  of  the  sP  phase  in  addition  to  the  pP  phase  in  the  test 
results  has  proved  useful  in  depth  determinations. 
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As  was  anticipated,  the  use  of  continental-size  arrays 
with  the  corresponding  increased  moveout  for  the  pP  phase  and 
the  large  spatial  separation  of  stations  does  serve  to  reduce  the 
coda-correlation  among  seismometers.  This  eliminates  the  presence 
of  confusing,  multiple  peaks  in  the  test  statistic,  which  were 
present  to  some  extent  when  the  automatic  test  was  performed  with 
LASA-Montana  data. 

For  two  events,  pre-filtering  of  the  data  was  found  to 
be  necessary.  After  filtering,  the  test  results  improved  vastly. 
The  erroneous  result  obtained  with  unfiltered  data  for  the 
19  July  1962  event  was  corrected  when  the  test  was  applied  to 
filtered  data,  yielding  clear  peaks  for  both  pF  and  sP.  Unfiltered 
data  for  the  6  October  1962  event  yielded  no  test  statistic  peaks, 
but  the  filtered  data  yielded  two  consistent  peaks  corresponding 
to  the  pP  and  sP  phase  arrivals. 

The  alignment  time  calculations  have  proved  adequate  for 
the  events  tested.  At  the  time  the  test  was  developed,  only  the 
JB  travel-time  tables  were  available  to  us.  It  was  therefore 
necessary  to  use  a  simple  model  to  calculate  sufficiently  precise 
P-pP  differential  travel  times.  More  recently,  we  have  obtained 
the  1966  Herrin  tables  [16],  which  include  pP  travel-time  tables 
that  would  be  adequate  for  the  automatic  pP  test.  Comparison  of 
differential  P-pP  arrival  times  obtained  from  our  model  and 
Herrin's  tables  indicate  that  the  discrepancies  are  probably  not 
important.  We  have  checked  two  ranges  for  each  of  two  depths 
corresponding  to  events  discussed  above.  In  one  case  both  our 
model  and  the  Herrin  table  predicts  the  same  differential  arrival 
time  for  a  depth  of  155  km.  In  the  other  case  the  differential 
arrival  time  predicted  by  the  model  for  a  depth  of  101  km  corres¬ 
ponds  to  a  depth  of  113  km  in  the  Herrin  tables.  Since  the  dis¬ 
crepancies  are  not  large,  a  conversion  of  the  test  to  take  ad¬ 
vantage  of  these  tables  does  not,  at  this  time,  seem  to  be  worth 
the  effort. 
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FIGURE  6.2  (  A-B) 

EXTENDED  AUTOMATIC  p  P  TEST 
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FIGURE  6.2  (  C * D  ) 
EXTENDED  AUTOMATIC  p  P  TEST 
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FIGURE  6.3  (  A-B) 

EXTENDED  AUTOMATIC  p  P  TEST 
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FIGURE  6.4 

SEISMOGRAMS  FOR  12  SEPTEMBER  1962  EVENT 


z 

Id 

> 

Id 

CNJ 

CD 

0> 


25 

0> 

Id  ^ 
00  £ 
3  O 
O  u- 

^  CD 


< 

OC 

o 

o 

2 

cn 

Id 

c/> 


-84 


85- 


SEISMOGRAMS  FOR  10 
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FIGURE  6.7 

SEISMOGRAMS  FOR  6  JANUARY  1964  EVENT 


APPENDIX  A 


VARIANCE  OF  INDIVIDUAL  ARRIVAL  TIME  ESTIMATES 


As  indicated  in  Section  III,  one  important  technique  for 
epicenter  estimation  begins  with  individual  arrival  time  estimates 
for  each  record.  The  variance  of  these  estimates  is  derived  in 
this  appendix  for  each  of  three  different  procedures.  In  each 
case  the  derived  results  apply  only  to  the  case  of  large  signal- 
to-noise  ratios  (and  hence  small  errors). 

For  the  purposes  of  this  derivation  the  assumptions  of 
Section  III  may  be  summarized  as  follows.  The  available  record, 
x(t) ,  is  given  by 

x(t)  =  s(t-x)  +  n(t)  (A-l) 


where  n(t)  is  stationary  and  has  power  density  spectrum  S  (f), 
and  s(t),  the  signal,  is  given  by 


and 


s(t)  =  sin2irf  t  -T/2  <  t  <  T/2 

=  0  1 1 1  >  T/2 

2lTf0  T  w0 


( A-2 ) 


(A-3) 


It  is  desired  to  estimate  ■?  from  x(t).  The  assumption  of  large 
signal-to-noise  ratio  implies  that  t  can  be  estimated  within  a 
fraction  of  T  by  a  cursory  examination  of  x(t).  Therefore,  the 
three  techniques  to  be  discussed  could  be  regarded  as  means  of 
refining  the  estimate  of  t  over  that  which  is  obvious  from  the 
unprocessed  record. 


A-l 


near 


ZERO  CROSSING  TIME  PICK 

In  this  case  t  is  defined  by  the  zero  crossing  of  x(t) 
Expanding  x(t)  about  t  we  have 

x(t )  =  x ( t )  +  ( t-i ) x ' ( t )+ .  .  .  (A- 4) 

=  s(0)+n(t)+(t-T)s' (0)+(t-i)n' (t)+. . . 
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Recognizing  that  s(0)  is  zero  and  defining  t  as  the  point  where 
x(t)  is  zero  yields 


0  =  n( t )+( t-t ) s ’ ( 0 )+( t-t )n ' ( t ) 
or 

i  -  T  n(x) 

T  "  T  "  s~(  0  )+n  ’T  t  ) 


(A-5) 
(A -6) 


Assuming  that  n'(-r)  may  be  neglected  compared  to  s'(0), 
yields 


this 

(A-7) 


Therefore,  the  variance  of  the  estimate  is  given  by 


v 


ar[ t ]  =  var[n(T>]  =  -i  j"  S  (f)df 
[s * ( o ) ;r  -« 


(A-8) 


A-2  MAXIMUM  VALUE  TIME  PICK 

Another  approach  to  estimating  the  arrival  time  would 
be  to  observe  the  time  near  t  +  T/4  where  a  maximum  in  x(t)  occurs. 
It  is  easily  seen  that  this  problem  is  equivalent  to  the  zero 
crossing  problem  just  considered  if  both  signal  and  noise  are 
replaced  by  their  derivatives.  We  have  therefore 


var[ t ] 


var[n 1 ( t ) ]  _  1 

[s"(T/iO]2  <3 


I 


(  2  irf )  2Sn  (  f )  df 


S  ( f ) df 
n 


(A-9) 


A-3  CROSS-CORRELATION  TIME  PICK 

The  third  method  of  estimating  arrival  times  is  more 
sensitive  to  variations  of  s(t)  from  its  assumed  shape.  In  the 
previous  two  methods,  only  the  behavior  of  s(t)  near  zero,  or 
near  T/^l  mattered,  and  therefore  the  analysis  might  still  be  valid 
even  if  s(t)  differed  in  several  respects  from  its  assumed  shape. 
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The  correlation  method  of  estimating  arrival  times  is  more  sensi¬ 
tive  to  variations  in  s(t),  and  might  deteriorate  badly  if  the 
received  signal  differed  appreciably  from  the  assumed  s(t). 

Basically  this  method  consists  of  choosing  t  as  the 
time  at  which  the  cross-correlation  function  of  x(t)  and  s(t) 
reaches  its  maximum.  This  is  equivalent  to  choosing  i  as  the 
time  of  maximum  output  for  a  linear  time  invariant  system  with 
x(t)  as  input  and  an  impulse  response  h(t)  given  by 

h(t)  -  s ( — t )  (A-10) 


This,  in  turn,  is  equivalent  to  looking  for  the  zero  crossing 
in  the  output  of  a  filter  with  impulse  response 

h'  (t)  =  -s ' (-t)  ( A-ll) 

Defining  y(t)  as  the  output  of  the  latter  system,  using  8  to 
indicate  convolution,  using  subscripts  s  and  n  to  indicate  the 
appropriate  components  in  y(t)  we  have 


y(t)  =  x(t)  8  h'(t)  =  -x(t)  0  s'(-t) 

=  -  s(t-x)  &  s'(-t)  -  n(t)  &  s 1 ( -t )  (A-12) 

V. _ _ _ J  V _ _ _ J 


=  ys(t-i)  +  yn(t) 

The  problem  is  now  in  the  form  of  the  zero  crossing  problem 
discussed  earlier.  That  is,  from  Equation  (A-7),  the  estimate 
is  given  by 


;  =  1 

It  only  remains  to  evaluate  y'(0)  and  var[y  (x)] 

S  II 

the  appropriate  substitution  in  Equation  (A-8). 


(A-13) 

before  making 


y  '  ( t )  =  +s ( t )  &  s" (-t ) 

o 

=  / s ( z ) s" ( - ( t-z ) ) dz 

y  '  (  0  )  =  /s(z)s"(z)dz 
s 


i  2,p  _ 

2W0 1  nu)0 


( A-l^ ) 
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The  power  density  spectrum  for  y  (t)  is  easily  obtained  as  the 
product  of  S  (f)  and  the  magnitude  squared  of  the  system  func¬ 
tion  corresponding  to  h'(t). 


S  (f)  =  S  (f)|/h'(t)e  j2lTftdt|2 
yn 

=  sn(f)  I /-s' (-t)e_,!2lTftdt|2 


f  f 

=  S  (f)[2sin^  •  -9--°-p-]2  (A-15) 

0 

The  detailed  calculations  required  to  obtain  the  last  line  of 
Equation  (A-15)  are  straightforward.  Finally,  we  have 


var[y  (  t  )  ]  =  fs  (f)df  =  fS  (f)[2sin^-  •  3— }-j]2df  (A-l6) 

yn  n  *0 

Substituting  Equations  ( A— 1 4 )  and  (A-l6)  into  Equation  (A-8), 
we  have 


var[y  (  t  )] 

var[  t  ]  =  - - — 3- 

Cy;(o)]2 


ff 


2  2 

f  -f 


°]2 


Sn(f)df 


(A-17) 
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APPENDIX  B 


EQUIVALENCE  OF  BEAMFORMING  AND  PLANE  WAVE  FIT 


B-l  CORRELATION  TIME  PICKS  IN  PRESENCE  OF  NOISE 

Referring  to  the  definitions  in  Section  III  of  the  text 
and  assuming  that  the  in  Equation  (3.1)  are  zero,  we  have  the 
following  expression  for  the  energy  in  the  output  of  the  array  as 
a  function  of  the  estimate  m. 


f(m)  =  f(  l  x .  ( t  +  im))2dt  (B-l) 

i  =  l  1 

N  N 

=11  l  [ s ( t-im+im) +n . (t+im) ][s(t-jm+jm) 

1=1  j  =  l 

+nj ( t+jm) ]dt 

It  is  assumed  that  the  limits  of  the  integral  roughly  bound 
the  duration  of  the  signal  terms;  this  will  be  possible  since 
a  large  signal-to-noise  ratio  is  assumed.  We  are  interested  in 
choosing  m  to  maximize  f (m) .  In  doing  so  we  shall  neglect  second- 
order  noise  contributions.  Discarding  these  terms  and  changing 
the  variables  of  integration 

N  N 

f (m)  =  l  £  js(t)s(t+(m-m) ( j-i) )dt 
1=1  j=l 

+ js ( t ) n . ( t+m( j -i ) +im)dt 

J 

+ Js( t)ni(t-m( j-i)+jm)dt  (B-2) 

In  order  to  find  the  m  for  which  f(m)  is  a  maximum,  we  solve 
for  f ' (m)  =0. 

f ' (m)  =  l  (j-i)  js(t)s’ (t+(m-m) ( j-i) )dt 

i j j  (B-3) 

+ Jdts ( t )  l  [ ( j-i)nl (t+m( j-i)+im)-( j-i)n! (t-m( j-i)+jm) ] 

i,j  J 


-91- 


Expanding  f'(m)  about  m, 

f'(m)  =  f'(m)  +  (m-m) f " (m)+ . . . 

Equating  this  to  zero  yields 

b,  _  „  l  -f ' (m) 
m  "  m  +  f"(m) 

Examining  Equation  ( B— 3 )  we  have 

f'(m)  =  0+/dts(t )  l  ( j-i ) [n' (t+jm)-n! (t+im) ] 

i  3  j  J 

•  n2 


(B-H) 


(B-5) 


( B— 6 ) 


4  \2 


f"(m)  =  l  (j-i)  /s(t)s"(t)dt  +  /dts(t)  £  (j-i) 
i  s  j  i ,  j 

•  [ n j  (t+jm)+n'.'(t+im)  ]  (B-7) 

In  Equation  (B-7),  which  gives  the  denominator  for  Equation  (B-5), 
we  neglect  the  noise  term  relative  to  the  signal  term.  Referring 
to  the  derivation  of  correlation  time  picks  (Equation  A-l^)  we 
may  recognize  the  signal  term 

' 2 'J*2 

0 


f "  (m)  =  [  l  ( j-i)2]/s(t)s"  (t)dt  =  y'(0)  1 — — £ - — 


=  —  TT  Cl) 


n2(n2-d 


0 


(B-8) 


In  Equation  (B-6)  it  is  useful  to  rewrite  the  double  sum  in¬ 
volving  the  noise  terms;  we  define  an  equivalent  noise  n^(t) 
by  this  summation: 

n^(t)  “  L  (j-i)  [nl  (t  +  jm)-nn!  (t  +  im)  ] 


i  j  j 

r 

L 

i  j  j 


/  r jnl -jn! -in !+in ! 1 
L  J  i  J  i 


=  2[N^in|  -  N(-f+1-)  l  n '  ] 
i  i 

N  M. 

=  2N[  l  (i  -  ^)ni] 

i  =  l  1 


(B-9) 
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We  are  now  ready  to  relate  the  beamforming  estimate  of  m 
to  that  which  would  be  obtained  by  Individual  correlation  time 
picks  followed  by  a  least  squares  slope  fit.  To  do  this,  recall 
Equation  ( A— 1 3 )  from  the  correlation  time  pick  discussion  which 
is  repeated  here  with  the  addition  of  a  subscript. 


Ti  Ti  " 


>'n1(Ti) 


( B-10 ) 


where  y  (t)  =  -n.(t)  8  s'(-t) 

1 

=  -/n1(t-z)s ' (-z)dz 
=-/n . ( t+z ' )s' ( z  T )  dz  ' 
=  +/s (z)n! (t+z)dz 

y  '  (  0  )  =  /s(z)s"(z)dz  =  -no >n 


and 


Assuming  the  true  arrival  times,  are  given  by 


yn  (ri)  =  +/s(z)n^(z+im)dz 


( B-ll ) 

( B-12 ) 

(B-13) 

(B-14) 


Substituting  Equations  (B-6),  (B-8)  and  ( B— 9 )  into  Equation  (B-5) 
we  have 

-/dts(t)n' (t) 


m  =  m  + 


N2(N"-1) 


y’(°) 


-2N£(i-  ^)/dts(t)n'(t+im) 


=  m  + 


=  m  + 


N2(N2-1) 

5 

y(1_  H±1)S 
L  2  K 

N(N2-1) 

12 


(B-15) 


(B-l6) 


y’(°) 


(B-17) 
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Equation  (B-17)  is  identical  to  Equation  (3-10)  of  Section  III 

1“  In 

if  yn./y*  is  identified  as  the  random  component  of  the  i 
1  s 

arrival-time  measurement.  And,  by  Equation  ( A— 1 3 )  of  Appendix  A, 
yn./y '  is  exactly  the  random  component  of  the  i^1  arrival  time 
measurement  when  the  correlation  method  of  arrival  time  esti¬ 
mation  is  used.  Therefore,  in  the  case  of  sufficiently  large 
signal-to-noise  ratios,  the  beamforming  method  and  the  corre¬ 
lation  time  pick  followed  by  a  least  squares  slope  fit  are 
identical.  That  is,  they  not  only  perform  equally  well,  they 
yield  the  same  answer. 

B-2  NOISELESS  CASE  WITH  RANDOM  TIMING  ERRORS 

In  the  noiseless  case  with  random  timing  errors,  we 

have 

x.(t)  =  s(t-im+T.)  ( B— 1 8 ) 

and  by  Equation  (3-10)  the  plane  wave  fit  yields 

v/.  N+lx 
L  ( 1  2  T  i 

m  =•  m  -  — - = - -  (B-19) 

N(N  -1) 

If,  as  in  the  previous  section,  we  consider 
f(m)  =  J  ( £x^  ( t  +  im))2dt 

-  f  l  s[ t-i ( m-m)+ r . ]s[ t-j (m-m) + t  . ]dt  (B-20) 

i  >  j  J 

and  choose  rh  to  maximize  the  resulting  expression,  we  obtain 

f'(m)  =  j  l  2is ' ( t-i (m-m)+ t . ) s ( t-j ( m-m) + t . ) dt 

i  j  j  1  J 

=  2  £  /is (t )s ' (t-( i-j ) (m-m)+( t . -r , ) )dt  (B-21) 

i,j  1  J 

Assuming  that  ( i- j ) ( m-m ) - ( t . - r . )  is  small  and  recognizing  that 

u  J 

if  it  were  zero  the  integral  would  be  zero,  we  have 

f'(rh)  -  2  l  i[  (i-j  )(m-m)-(i.-T  ,)]/s(t)s"(t.)dt  =  0  (B-22) 

ij  1  J 

and  therefore. 
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( B— 2  3 ) 


(m-m)  l  (i2-ij)  =  £  i(t.-r.) 

ij  i,J  J 


which  reduces  to 


m  =  m 


ICi  - 

i 


N+lx 
— )xi 


N(N- -1) 

12 


in  agreement  with  Equation  (B-19). 


(B-24) 
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